home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgesvd.f < prev    next >
Text File  |  1996-07-19  |  134KB  |  3,409 lines
  1.       SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
  2.      $                   WORK, LWORK, INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOBU, JOBVT
  11.       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
  15.      $                   VT( LDVT, * ), WORK( * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  DGESVD computes the singular value decomposition (SVD) of a real
  22. *  M-by-N matrix A, optionally computing the left and/or right singular
  23. *  vectors. The SVD is written
  24. *
  25. *       A = U * SIGMA * transpose(V)
  26. *
  27. *  where SIGMA is an M-by-N matrix which is zero except for its
  28. *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  29. *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  30. *  are the singular values of A; they are real and non-negative, and
  31. *  are returned in descending order.  The first min(m,n) columns of
  32. *  U and V are the left and right singular vectors of A.
  33. *
  34. *  Note that the routine returns V**T, not V.
  35. *
  36. *  Arguments
  37. *  =========
  38. *
  39. *  JOBU    (input) CHARACTER*1
  40. *          Specifies options for computing all or part of the matrix U:
  41. *          = 'A':  all M columns of U are returned in array U:
  42. *          = 'S':  the first min(m,n) columns of U (the left singular
  43. *                  vectors) are returned in the array U;
  44. *          = 'O':  the first min(m,n) columns of U (the left singular
  45. *                  vectors) are overwritten on the array A;
  46. *          = 'N':  no columns of U (no left singular vectors) are
  47. *                  computed.
  48. *
  49. *  JOBVT   (input) CHARACTER*1
  50. *          Specifies options for computing all or part of the matrix
  51. *          V**T:
  52. *          = 'A':  all N rows of V**T are returned in the array VT;
  53. *          = 'S':  the first min(m,n) rows of V**T (the right singular
  54. *                  vectors) are returned in the array VT;
  55. *          = 'O':  the first min(m,n) rows of V**T (the right singular
  56. *                  vectors) are overwritten on the array A;
  57. *          = 'N':  no rows of V**T (no right singular vectors) are
  58. *                  computed.
  59. *
  60. *          JOBVT and JOBU cannot both be 'O'.
  61. *
  62. *  M       (input) INTEGER
  63. *          The number of rows of the input matrix A.  M >= 0.
  64. *
  65. *  N       (input) INTEGER
  66. *          The number of columns of the input matrix A.  N >= 0.
  67. *
  68. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  69. *          On entry, the M-by-N matrix A.
  70. *          On exit,
  71. *          if JOBU = 'O',  A is overwritten with the first min(m,n)
  72. *                          columns of U (the left singular vectors,
  73. *                          stored columnwise);
  74. *          if JOBVT = 'O', A is overwritten with the first min(m,n)
  75. *                          rows of V**T (the right singular vectors,
  76. *                          stored rowwise);
  77. *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
  78. *                          are destroyed.
  79. *
  80. *  LDA     (input) INTEGER
  81. *          The leading dimension of the array A.  LDA >= max(1,M).
  82. *
  83. *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
  84. *          The singular values of A, sorted so that S(i) >= S(i+1).
  85. *
  86. *  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
  87. *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
  88. *          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
  89. *          if JOBU = 'S', U contains the first min(m,n) columns of U
  90. *          (the left singular vectors, stored columnwise);
  91. *          if JOBU = 'N' or 'O', U is not referenced.
  92. *
  93. *  LDU     (input) INTEGER
  94. *          The leading dimension of the array U.  LDU >= 1; if
  95. *          JOBU = 'S' or 'A', LDU >= M.
  96. *
  97. *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
  98. *          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
  99. *          V**T;
  100. *          if JOBVT = 'S', VT contains the first min(m,n) rows of
  101. *          V**T (the right singular vectors, stored rowwise);
  102. *          if JOBVT = 'N' or 'O', VT is not referenced.
  103. *
  104. *  LDVT    (input) INTEGER
  105. *          The leading dimension of the array VT.  LDVT >= 1; if
  106. *          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
  107. *
  108. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  109. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
  110. *          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
  111. *          superdiagonal elements of an upper bidiagonal matrix B
  112. *          whose diagonal is in S (not necessarily sorted). B
  113. *          satisfies A = U * B * VT, so it has the same singular values
  114. *          as A, and singular vectors related by U and VT.
  115. *
  116. *  LWORK   (input) INTEGER
  117. *          The dimension of the array WORK. LWORK >= 1.
  118. *          LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)-4).
  119. *          For good performance, LWORK should generally be larger.
  120. *
  121. *  INFO    (output) INTEGER
  122. *          = 0:  successful exit.
  123. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  124. *          > 0:  if DBDSQR did not converge, INFO specifies how many
  125. *                superdiagonals of an intermediate bidiagonal form B
  126. *                did not converge to zero. See the description of WORK
  127. *                above for details.
  128. *
  129. *  =====================================================================
  130. *
  131. *     .. Parameters ..
  132.       DOUBLE PRECISION   ZERO, ONE
  133.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  134. *     ..
  135. *     .. Local Scalars ..
  136.       LOGICAL            WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS, WNTVA,
  137.      $                   WNTVAS, WNTVN, WNTVO, WNTVS
  138.       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IR, ISCL,
  139.      $                   ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
  140.      $                   MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
  141.      $                   NRVT, WRKBL
  142.       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
  143. *     ..
  144. *     .. Local Arrays ..
  145.       DOUBLE PRECISION   DUM( 1 )
  146. *     ..
  147. *     .. External Subroutines ..
  148.       EXTERNAL           DBDSQR, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
  149.      $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
  150.      $                   XERBLA
  151. *     ..
  152. *     .. External Functions ..
  153.       LOGICAL            LSAME
  154.       INTEGER            ILAENV
  155.       DOUBLE PRECISION   DLAMCH, DLANGE
  156.       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
  157. *     ..
  158. *     .. Intrinsic Functions ..
  159.       INTRINSIC          MAX, MIN, SQRT
  160. *     ..
  161. *     .. Executable Statements ..
  162. *
  163. *     Test the input arguments
  164. *
  165.       INFO = 0
  166.       MINMN = MIN( M, N )
  167.       MNTHR = ILAENV( 6, 'DGESVD', JOBU // JOBVT, M, N, 0, 0 )
  168.       WNTUA = LSAME( JOBU, 'A' )
  169.       WNTUS = LSAME( JOBU, 'S' )
  170.       WNTUAS = WNTUA .OR. WNTUS
  171.       WNTUO = LSAME( JOBU, 'O' )
  172.       WNTUN = LSAME( JOBU, 'N' )
  173.       WNTVA = LSAME( JOBVT, 'A' )
  174.       WNTVS = LSAME( JOBVT, 'S' )
  175.       WNTVAS = WNTVA .OR. WNTVS
  176.       WNTVO = LSAME( JOBVT, 'O' )
  177.       WNTVN = LSAME( JOBVT, 'N' )
  178.       MINWRK = 1
  179. *
  180.       IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
  181.          INFO = -1
  182.       ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
  183.      $         ( WNTVO .AND. WNTUO ) ) THEN
  184.          INFO = -2
  185.       ELSE IF( M.LT.0 ) THEN
  186.          INFO = -3
  187.       ELSE IF( N.LT.0 ) THEN
  188.          INFO = -4
  189.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  190.          INFO = -6
  191.       ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
  192.          INFO = -9
  193.       ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
  194.      $         ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
  195.          INFO = -11
  196.       END IF
  197. *
  198. *     Compute workspace
  199. *      (Note: Comments in the code beginning "Workspace:" describe the
  200. *       minimal amount of workspace needed at that point in the code,
  201. *       as well as the preferred amount for good performance.
  202. *       NB refers to the optimal block size for the immediately
  203. *       following subroutine, as returned by ILAENV.)
  204. *
  205.       IF( INFO.EQ.0 .AND. LWORK.GE.1 .AND. M.GT.0 .AND. N.GT.0 ) THEN
  206.          IF( M.GE.N ) THEN
  207. *
  208. *           Compute space needed for DBDSQR
  209. *
  210.             BDSPAC = MAX( 3*N, 5*N-4 )
  211.             IF( M.GE.MNTHR ) THEN
  212.                IF( WNTUN ) THEN
  213. *
  214. *                 Path 1 (M much larger than N, JOBU='N')
  215. *
  216.                   MAXWRK = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
  217.      $                     -1 )
  218.                   MAXWRK = MAX( MAXWRK, 3*N+2*N*
  219.      $                     ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  220.                   IF( WNTVO .OR. WNTVAS )
  221.      $               MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
  222.      $                        ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  223.                   MAXWRK = MAX( MAXWRK, BDSPAC )
  224.                   MINWRK = MAX( 4*N, BDSPAC )
  225.                   MAXWRK = MAX( MAXWRK, MINWRK )
  226.                ELSE IF( WNTUO .AND. WNTVN ) THEN
  227. *
  228. *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
  229. *
  230.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  231.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
  232.      $                    N, N, -1 ) )
  233.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  234.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  235.                   WRKBL = MAX( WRKBL, 3*N+N*
  236.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  237.                   WRKBL = MAX( WRKBL, BDSPAC )
  238.                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
  239.                   MINWRK = MAX( 3*N+M, BDSPAC )
  240.                   MAXWRK = MAX( MAXWRK, MINWRK )
  241.                ELSE IF( WNTUO .AND. WNTVAS ) THEN
  242. *
  243. *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
  244. *                 'A')
  245. *
  246.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  247.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
  248.      $                    N, N, -1 ) )
  249.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  250.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  251.                   WRKBL = MAX( WRKBL, 3*N+N*
  252.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  253.                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
  254.      $                    ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  255.                   WRKBL = MAX( WRKBL, BDSPAC )
  256.                   MAXWRK = MAX( N*N+WRKBL, N*N+M*N+N )
  257.                   MINWRK = MAX( 3*N+M, BDSPAC )
  258.                   MAXWRK = MAX( MAXWRK, MINWRK )
  259.                ELSE IF( WNTUS .AND. WNTVN ) THEN
  260. *
  261. *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
  262. *
  263.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  264.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
  265.      $                    N, N, -1 ) )
  266.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  267.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  268.                   WRKBL = MAX( WRKBL, 3*N+N*
  269.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  270.                   WRKBL = MAX( WRKBL, BDSPAC )
  271.                   MAXWRK = N*N + WRKBL
  272.                   MINWRK = MAX( 3*N+M, BDSPAC )
  273.                   MAXWRK = MAX( MAXWRK, MINWRK )
  274.                ELSE IF( WNTUS .AND. WNTVO ) THEN
  275. *
  276. *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
  277. *
  278.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  279.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
  280.      $                    N, N, -1 ) )
  281.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  282.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  283.                   WRKBL = MAX( WRKBL, 3*N+N*
  284.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  285.                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
  286.      $                    ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  287.                   WRKBL = MAX( WRKBL, BDSPAC )
  288.                   MAXWRK = 2*N*N + WRKBL
  289.                   MINWRK = MAX( 3*N+M, BDSPAC )
  290.                   MAXWRK = MAX( MAXWRK, MINWRK )
  291.                ELSE IF( WNTUS .AND. WNTVAS ) THEN
  292. *
  293. *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
  294. *                 'A')
  295. *
  296.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  297.                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
  298.      $                    N, N, -1 ) )
  299.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  300.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  301.                   WRKBL = MAX( WRKBL, 3*N+N*
  302.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  303.                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
  304.      $                    ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  305.                   WRKBL = MAX( WRKBL, BDSPAC )
  306.                   MAXWRK = N*N + WRKBL
  307.                   MINWRK = MAX( 3*N+M, BDSPAC )
  308.                   MAXWRK = MAX( MAXWRK, MINWRK )
  309.                ELSE IF( WNTUA .AND. WNTVN ) THEN
  310. *
  311. *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
  312. *
  313.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  314.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
  315.      $                    M, N, -1 ) )
  316.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  317.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  318.                   WRKBL = MAX( WRKBL, 3*N+N*
  319.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  320.                   WRKBL = MAX( WRKBL, BDSPAC )
  321.                   MAXWRK = N*N + WRKBL
  322.                   MINWRK = MAX( 3*N+M, BDSPAC )
  323.                   MAXWRK = MAX( MAXWRK, MINWRK )
  324.                ELSE IF( WNTUA .AND. WNTVO ) THEN
  325. *
  326. *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
  327. *
  328.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  329.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
  330.      $                    M, N, -1 ) )
  331.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  332.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  333.                   WRKBL = MAX( WRKBL, 3*N+N*
  334.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  335.                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
  336.      $                    ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  337.                   WRKBL = MAX( WRKBL, BDSPAC )
  338.                   MAXWRK = 2*N*N + WRKBL
  339.                   MINWRK = MAX( 3*N+M, BDSPAC )
  340.                   MAXWRK = MAX( MAXWRK, MINWRK )
  341.                ELSE IF( WNTUA .AND. WNTVAS ) THEN
  342. *
  343. *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
  344. *                 'A')
  345. *
  346.                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
  347.                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
  348.      $                    M, N, -1 ) )
  349.                   WRKBL = MAX( WRKBL, 3*N+2*N*
  350.      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
  351.                   WRKBL = MAX( WRKBL, 3*N+N*
  352.      $                    ILAENV( 1, 'DORGBR', 'Q', N, N, N, -1 ) )
  353.                   WRKBL = MAX( WRKBL, 3*N+( N-1 )*
  354.      $                    ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  355.                   WRKBL = MAX( WRKBL, BDSPAC )
  356.                   MAXWRK = N*N + WRKBL
  357.                   MINWRK = MAX( 3*N+M, BDSPAC )
  358.                   MAXWRK = MAX( MAXWRK, MINWRK )
  359.                END IF
  360.             ELSE
  361. *
  362. *              Path 10 (M at least N, but not much larger)
  363. *
  364.                MAXWRK = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N,
  365.      $                  -1, -1 )
  366.                IF( WNTUS .OR. WNTUO )
  367.      $            MAXWRK = MAX( MAXWRK, 3*N+N*
  368.      $                     ILAENV( 1, 'DORGBR', 'Q', M, N, N, -1 ) )
  369.                IF( WNTUA )
  370.      $            MAXWRK = MAX( MAXWRK, 3*N+M*
  371.      $                     ILAENV( 1, 'DORGBR', 'Q', M, M, N, -1 ) )
  372.                IF( .NOT.WNTVN )
  373.      $            MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
  374.      $                     ILAENV( 1, 'DORGBR', 'P', N, N, N, -1 ) )
  375.                MAXWRK = MAX( MAXWRK, BDSPAC )
  376.                MINWRK = MAX( 3*N+M, BDSPAC )
  377.                MAXWRK = MAX( MAXWRK, MINWRK )
  378.             END IF
  379.          ELSE
  380. *
  381. *           Compute space needed for DBDSQR
  382. *
  383.             BDSPAC = MAX( 3*M, 5*M-4 )
  384.             IF( N.GE.MNTHR ) THEN
  385.                IF( WNTVN ) THEN
  386. *
  387. *                 Path 1t(N much larger than M, JOBVT='N')
  388. *
  389.                   MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
  390.      $                     -1 )
  391.                   MAXWRK = MAX( MAXWRK, 3*M+2*M*
  392.      $                     ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  393.                   IF( WNTUO .OR. WNTUAS )
  394.      $               MAXWRK = MAX( MAXWRK, 3*M+M*
  395.      $                        ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  396.                   MAXWRK = MAX( MAXWRK, BDSPAC )
  397.                   MINWRK = MAX( 4*M, BDSPAC )
  398.                   MAXWRK = MAX( MAXWRK, MINWRK )
  399.                ELSE IF( WNTVO .AND. WNTUN ) THEN
  400. *
  401. *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
  402. *
  403.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  404.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
  405.      $                    N, M, -1 ) )
  406.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  407.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  408.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  409.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  410.                   WRKBL = MAX( WRKBL, BDSPAC )
  411.                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
  412.                   MINWRK = MAX( 3*M+N, BDSPAC )
  413.                   MAXWRK = MAX( MAXWRK, MINWRK )
  414.                ELSE IF( WNTVO .AND. WNTUAS ) THEN
  415. *
  416. *                 Path 3t(N much larger than M, JOBU='S' or 'A',
  417. *                 JOBVT='O')
  418. *
  419.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  420.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
  421.      $                    N, M, -1 ) )
  422.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  423.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  424.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  425.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  426.                   WRKBL = MAX( WRKBL, 3*M+M*
  427.      $                    ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  428.                   WRKBL = MAX( WRKBL, BDSPAC )
  429.                   MAXWRK = MAX( M*M+WRKBL, M*M+M*N+M )
  430.                   MINWRK = MAX( 3*M+N, BDSPAC )
  431.                   MAXWRK = MAX( MAXWRK, MINWRK )
  432.                ELSE IF( WNTVS .AND. WNTUN ) THEN
  433. *
  434. *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
  435. *
  436.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  437.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
  438.      $                    N, M, -1 ) )
  439.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  440.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  441.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  442.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  443.                   WRKBL = MAX( WRKBL, BDSPAC )
  444.                   MAXWRK = M*M + WRKBL
  445.                   MINWRK = MAX( 3*M+N, BDSPAC )
  446.                   MAXWRK = MAX( MAXWRK, MINWRK )
  447.                ELSE IF( WNTVS .AND. WNTUO ) THEN
  448. *
  449. *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
  450. *
  451.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  452.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
  453.      $                    N, M, -1 ) )
  454.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  455.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  456.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  457.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  458.                   WRKBL = MAX( WRKBL, 3*M+M*
  459.      $                    ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  460.                   WRKBL = MAX( WRKBL, BDSPAC )
  461.                   MAXWRK = 2*M*M + WRKBL
  462.                   MINWRK = MAX( 3*M+N, BDSPAC )
  463.                   MAXWRK = MAX( MAXWRK, MINWRK )
  464.                ELSE IF( WNTVS .AND. WNTUAS ) THEN
  465. *
  466. *                 Path 6t(N much larger than M, JOBU='S' or 'A',
  467. *                 JOBVT='S')
  468. *
  469.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  470.                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
  471.      $                    N, M, -1 ) )
  472.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  473.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  474.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  475.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  476.                   WRKBL = MAX( WRKBL, 3*M+M*
  477.      $                    ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  478.                   WRKBL = MAX( WRKBL, BDSPAC )
  479.                   MAXWRK = M*M + WRKBL
  480.                   MINWRK = MAX( 3*M+N, BDSPAC )
  481.                   MAXWRK = MAX( MAXWRK, MINWRK )
  482.                ELSE IF( WNTVA .AND. WNTUN ) THEN
  483. *
  484. *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
  485. *
  486.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  487.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
  488.      $                    N, M, -1 ) )
  489.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  490.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  491.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  492.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  493.                   WRKBL = MAX( WRKBL, BDSPAC )
  494.                   MAXWRK = M*M + WRKBL
  495.                   MINWRK = MAX( 3*M+N, BDSPAC )
  496.                   MAXWRK = MAX( MAXWRK, MINWRK )
  497.                ELSE IF( WNTVA .AND. WNTUO ) THEN
  498. *
  499. *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
  500. *
  501.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  502.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
  503.      $                    N, M, -1 ) )
  504.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  505.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  506.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  507.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  508.                   WRKBL = MAX( WRKBL, 3*M+M*
  509.      $                    ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  510.                   WRKBL = MAX( WRKBL, BDSPAC )
  511.                   MAXWRK = 2*M*M + WRKBL
  512.                   MINWRK = MAX( 3*M+N, BDSPAC )
  513.                   MAXWRK = MAX( MAXWRK, MINWRK )
  514.                ELSE IF( WNTVA .AND. WNTUAS ) THEN
  515. *
  516. *                 Path 9t(N much larger than M, JOBU='S' or 'A',
  517. *                 JOBVT='A')
  518. *
  519.                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
  520.                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
  521.      $                    N, M, -1 ) )
  522.                   WRKBL = MAX( WRKBL, 3*M+2*M*
  523.      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
  524.                   WRKBL = MAX( WRKBL, 3*M+( M-1 )*
  525.      $                    ILAENV( 1, 'DORGBR', 'P', M, M, M, -1 ) )
  526.                   WRKBL = MAX( WRKBL, 3*M+M*
  527.      $                    ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  528.                   WRKBL = MAX( WRKBL, BDSPAC )
  529.                   MAXWRK = M*M + WRKBL
  530.                   MINWRK = MAX( 3*M+N, BDSPAC )
  531.                   MAXWRK = MAX( MAXWRK, MINWRK )
  532.                END IF
  533.             ELSE
  534. *
  535. *              Path 10t(N greater than M, but not much larger)
  536. *
  537.                MAXWRK = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N,
  538.      $                  -1, -1 )
  539.                IF( WNTVS .OR. WNTVO )
  540.      $            MAXWRK = MAX( MAXWRK, 3*M+M*
  541.      $                     ILAENV( 1, 'DORGBR', 'P', M, N, M, -1 ) )
  542.                IF( WNTVA )
  543.      $            MAXWRK = MAX( MAXWRK, 3*M+N*
  544.      $                     ILAENV( 1, 'DORGBR', 'P', N, N, M, -1 ) )
  545.                IF( .NOT.WNTUN )
  546.      $            MAXWRK = MAX( MAXWRK, 3*M+( M-1 )*
  547.      $                     ILAENV( 1, 'DORGBR', 'Q', M, M, M, -1 ) )
  548.                MAXWRK = MAX( MAXWRK, BDSPAC )
  549.                MINWRK = MAX( 3*M+N, BDSPAC )
  550.                MAXWRK = MAX( MAXWRK, MINWRK )
  551.             END IF
  552.          END IF
  553.          WORK( 1 ) = MAXWRK
  554.       END IF
  555. *
  556.       IF( LWORK.LT.MINWRK ) THEN
  557.          INFO = -13
  558.       END IF
  559.       IF( INFO.NE.0 ) THEN
  560.          CALL XERBLA( 'DGESVD', -INFO )
  561.          RETURN
  562.       END IF
  563. *
  564. *     Quick return if possible
  565. *
  566.       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  567.          IF( LWORK.GE.1 )
  568.      $      WORK( 1 ) = ONE
  569.          RETURN
  570.       END IF
  571. *
  572. *     Get machine constants
  573. *
  574.       EPS = DLAMCH( 'P' )
  575.       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
  576.       BIGNUM = ONE / SMLNUM
  577. *
  578. *     Scale A if max element outside range [SMLNUM,BIGNUM]
  579. *
  580.       ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
  581.       ISCL = 0
  582.       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  583.          ISCL = 1
  584.          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
  585.       ELSE IF( ANRM.GT.BIGNUM ) THEN
  586.          ISCL = 1
  587.          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
  588.       END IF
  589. *
  590.       IF( M.GE.N ) THEN
  591. *
  592. *        A has at least as many rows as columns. If A has sufficiently
  593. *        more rows than columns, first reduce using the QR
  594. *        decomposition (if sufficient workspace available)
  595. *
  596.          IF( M.GE.MNTHR ) THEN
  597. *
  598.             IF( WNTUN ) THEN
  599. *
  600. *              Path 1 (M much larger than N, JOBU='N')
  601. *              No left singular vectors to be computed
  602. *
  603.                ITAU = 1
  604.                IWORK = ITAU + N
  605. *
  606. *              Compute A=Q*R
  607. *              (Workspace: need 2*N, prefer N+N*NB)
  608. *
  609.                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  610.      $                      LWORK-IWORK+1, IERR )
  611. *
  612. *              Zero out below R
  613. *
  614.                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
  615.                IE = 1
  616.                ITAUQ = IE + N
  617.                ITAUP = ITAUQ + N
  618.                IWORK = ITAUP + N
  619. *
  620. *              Bidiagonalize R in A
  621. *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
  622. *
  623.                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  624.      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  625.      $                      IERR )
  626.                NCVT = 0
  627.                IF( WNTVO .OR. WNTVAS ) THEN
  628. *
  629. *                 If right singular vectors desired, generate P'.
  630. *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  631. *
  632.                   CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  633.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  634.                   NCVT = N
  635.                END IF
  636.                IWORK = IE + N
  637. *
  638. *              Perform bidiagonal QR iteration, computing right
  639. *              singular vectors of A in A if desired
  640. *              (Workspace: need BDSPAC)
  641. *
  642.                CALL DBDSQR( 'U', N, NCVT, 0, 0, S, WORK( IE ), A, LDA,
  643.      $                      DUM, 1, DUM, 1, WORK( IWORK ), INFO )
  644. *
  645. *              If right singular vectors desired in VT, copy them there
  646. *
  647.                IF( WNTVAS )
  648.      $            CALL DLACPY( 'F', N, N, A, LDA, VT, LDVT )
  649. *
  650.             ELSE IF( WNTUO .AND. WNTVN ) THEN
  651. *
  652. *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
  653. *              N left singular vectors to be overwritten on A and
  654. *              no right singular vectors to be computed
  655. *
  656.                IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
  657. *
  658. *                 Sufficient workspace for a fast algorithm
  659. *
  660.                   IR = 1
  661.                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
  662. *
  663. *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
  664. *
  665.                      LDWRKU = LDA
  666.                      LDWRKR = LDA
  667.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
  668. *
  669. *                    WORK(IU) is LDA by N, WORK(IR) is N by N
  670. *
  671.                      LDWRKU = LDA
  672.                      LDWRKR = N
  673.                   ELSE
  674. *
  675. *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
  676. *
  677.                      LDWRKU = ( LWORK-N*N-N ) / N
  678.                      LDWRKR = N
  679.                   END IF
  680.                   ITAU = IR + LDWRKR*N
  681.                   IWORK = ITAU + N
  682. *
  683. *                 Compute A=Q*R
  684. *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  685. *
  686.                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  687.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  688. *
  689. *                 Copy R to WORK(IR) and zero out below it
  690. *
  691.                   CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
  692.                   CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
  693.      $                         LDWRKR )
  694. *
  695. *                 Generate Q in A
  696. *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  697. *
  698.                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  699.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  700.                   IE = ITAU
  701.                   ITAUQ = IE + N
  702.                   ITAUP = ITAUQ + N
  703.                   IWORK = ITAUP + N
  704. *
  705. *                 Bidiagonalize R in WORK(IR)
  706. *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  707. *
  708.                   CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
  709.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  710.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  711. *
  712. *                 Generate left vectors bidiagonalizing R
  713. *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  714. *
  715.                   CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  716.      $                         WORK( ITAUQ ), WORK( IWORK ),
  717.      $                         LWORK-IWORK+1, IERR )
  718.                   IWORK = IE + N
  719. *
  720. *                 Perform bidiagonal QR iteration, computing left
  721. *                 singular vectors of R in WORK(IR)
  722. *                 (Workspace: need N*N+BDSPAC)
  723. *
  724.                   CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM, 1,
  725.      $                         WORK( IR ), LDWRKR, DUM, 1,
  726.      $                         WORK( IWORK ), INFO )
  727.                   IU = IE + N
  728. *
  729. *                 Multiply Q in A by left singular vectors of R in
  730. *                 WORK(IR), storing result in WORK(IU) and copying to A
  731. *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
  732. *
  733.                   DO 10 I = 1, M, LDWRKU
  734.                      CHUNK = MIN( M-I+1, LDWRKU )
  735.                      CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
  736.      $                           LDA, WORK( IR ), LDWRKR, ZERO,
  737.      $                           WORK( IU ), LDWRKU )
  738.                      CALL DLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
  739.      $                            A( I, 1 ), LDA )
  740.    10             CONTINUE
  741. *
  742.                ELSE
  743. *
  744. *                 Insufficient workspace for a fast algorithm
  745. *
  746.                   IE = 1
  747.                   ITAUQ = IE + N
  748.                   ITAUP = ITAUQ + N
  749.                   IWORK = ITAUP + N
  750. *
  751. *                 Bidiagonalize A
  752. *                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
  753. *
  754.                   CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
  755.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  756.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  757. *
  758. *                 Generate left vectors bidiagonalizing A
  759. *                 (Workspace: need 4*N, prefer 3*N+N*NB)
  760. *
  761.                   CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
  762.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  763.                   IWORK = IE + N
  764. *
  765. *                 Perform bidiagonal QR iteration, computing left
  766. *                 singular vectors of A in A
  767. *                 (Workspace: need BDSPAC)
  768. *
  769.                   CALL DBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM, 1,
  770.      $                         A, LDA, DUM, 1, WORK( IWORK ), INFO )
  771. *
  772.                END IF
  773. *
  774.             ELSE IF( WNTUO .AND. WNTVAS ) THEN
  775. *
  776. *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
  777. *              N left singular vectors to be overwritten on A and
  778. *              N right singular vectors to be computed in VT
  779. *
  780.                IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
  781. *
  782. *                 Sufficient workspace for a fast algorithm
  783. *
  784.                   IR = 1
  785.                   IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+LDA*N ) THEN
  786. *
  787. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
  788. *
  789.                      LDWRKU = LDA
  790.                      LDWRKR = LDA
  791.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+N )+N*N ) THEN
  792. *
  793. *                    WORK(IU) is LDA by N and WORK(IR) is N by N
  794. *
  795.                      LDWRKU = LDA
  796.                      LDWRKR = N
  797.                   ELSE
  798. *
  799. *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
  800. *
  801.                      LDWRKU = ( LWORK-N*N-N ) / N
  802.                      LDWRKR = N
  803.                   END IF
  804.                   ITAU = IR + LDWRKR*N
  805.                   IWORK = ITAU + N
  806. *
  807. *                 Compute A=Q*R
  808. *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  809. *
  810.                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  811.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  812. *
  813. *                 Copy R to VT, zeroing out below it
  814. *
  815.                   CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
  816.                   CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, VT( 2, 1 ),
  817.      $                         LDVT )
  818. *
  819. *                 Generate Q in A
  820. *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  821. *
  822.                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  823.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  824.                   IE = ITAU
  825.                   ITAUQ = IE + N
  826.                   ITAUP = ITAUQ + N
  827.                   IWORK = ITAUP + N
  828. *
  829. *                 Bidiagonalize R in VT, copying result to WORK(IR)
  830. *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  831. *
  832.                   CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
  833.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  834.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  835.                   CALL DLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
  836. *
  837. *                 Generate left vectors bidiagonalizing R in WORK(IR)
  838. *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  839. *
  840.                   CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  841.      $                         WORK( ITAUQ ), WORK( IWORK ),
  842.      $                         LWORK-IWORK+1, IERR )
  843. *
  844. *                 Generate right vectors bidiagonalizing R in VT
  845. *                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
  846. *
  847.                   CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  848.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  849.                   IWORK = IE + N
  850. *
  851. *                 Perform bidiagonal QR iteration, computing left
  852. *                 singular vectors of R in WORK(IR) and computing right
  853. *                 singular vectors of R in VT
  854. *                 (Workspace: need N*N+BDSPAC)
  855. *
  856.                   CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT, LDVT,
  857.      $                         WORK( IR ), LDWRKR, DUM, 1,
  858.      $                         WORK( IWORK ), INFO )
  859.                   IU = IE + N
  860. *
  861. *                 Multiply Q in A by left singular vectors of R in
  862. *                 WORK(IR), storing result in WORK(IU) and copying to A
  863. *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
  864. *
  865.                   DO 20 I = 1, M, LDWRKU
  866.                      CHUNK = MIN( M-I+1, LDWRKU )
  867.                      CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
  868.      $                           LDA, WORK( IR ), LDWRKR, ZERO,
  869.      $                           WORK( IU ), LDWRKU )
  870.                      CALL DLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
  871.      $                            A( I, 1 ), LDA )
  872.    20             CONTINUE
  873. *
  874.                ELSE
  875. *
  876. *                 Insufficient workspace for a fast algorithm
  877. *
  878.                   ITAU = 1
  879.                   IWORK = ITAU + N
  880. *
  881. *                 Compute A=Q*R
  882. *                 (Workspace: need 2*N, prefer N+N*NB)
  883. *
  884.                   CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  885.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  886. *
  887. *                 Copy R to VT, zeroing out below it
  888. *
  889.                   CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
  890.                   CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, VT( 2, 1 ),
  891.      $                         LDVT )
  892. *
  893. *                 Generate Q in A
  894. *                 (Workspace: need 2*N, prefer N+N*NB)
  895. *
  896.                   CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  897.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  898.                   IE = ITAU
  899.                   ITAUQ = IE + N
  900.                   ITAUP = ITAUQ + N
  901.                   IWORK = ITAUP + N
  902. *
  903. *                 Bidiagonalize R in VT
  904. *                 (Workspace: need 4*N, prefer 3*N+2*N*NB)
  905. *
  906.                   CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
  907.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  908.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  909. *
  910. *                 Multiply Q in A by left vectors bidiagonalizing R
  911. *                 (Workspace: need 3*N+M, prefer 3*N+M*NB)
  912. *
  913.                   CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  914.      $                         WORK( ITAUQ ), A, LDA, WORK( IWORK ),
  915.      $                         LWORK-IWORK+1, IERR )
  916. *
  917. *                 Generate right vectors bidiagonalizing R in VT
  918. *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  919. *
  920.                   CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  921.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  922.                   IWORK = IE + N
  923. *
  924. *                 Perform bidiagonal QR iteration, computing left
  925. *                 singular vectors of A in A and computing right
  926. *                 singular vectors of A in VT
  927. *                 (Workspace: need BDSPAC)
  928. *
  929.                   CALL DBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT, LDVT,
  930.      $                         A, LDA, DUM, 1, WORK( IWORK ), INFO )
  931. *
  932.                END IF
  933. *
  934.             ELSE IF( WNTUS ) THEN
  935. *
  936.                IF( WNTVN ) THEN
  937. *
  938. *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
  939. *                 N left singular vectors to be computed in U and
  940. *                 no right singular vectors to be computed
  941. *
  942.                   IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
  943. *
  944. *                    Sufficient workspace for a fast algorithm
  945. *
  946.                      IR = 1
  947.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  948. *
  949. *                       WORK(IR) is LDA by N
  950. *
  951.                         LDWRKR = LDA
  952.                      ELSE
  953. *
  954. *                       WORK(IR) is N by N
  955. *
  956.                         LDWRKR = N
  957.                      END IF
  958.                      ITAU = IR + LDWRKR*N
  959.                      IWORK = ITAU + N
  960. *
  961. *                    Compute A=Q*R
  962. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  963. *
  964.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  965.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  966. *
  967. *                    Copy R to WORK(IR), zeroing out below it
  968. *
  969.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ),
  970.      $                            LDWRKR )
  971.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  972.      $                            WORK( IR+1 ), LDWRKR )
  973. *
  974. *                    Generate Q in A
  975. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  976. *
  977.                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  978.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  979.                      IE = ITAU
  980.                      ITAUQ = IE + N
  981.                      ITAUP = ITAUQ + N
  982.                      IWORK = ITAUP + N
  983. *
  984. *                    Bidiagonalize R in WORK(IR)
  985. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  986. *
  987.                      CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
  988.      $                            WORK( IE ), WORK( ITAUQ ),
  989.      $                            WORK( ITAUP ), WORK( IWORK ),
  990.      $                            LWORK-IWORK+1, IERR )
  991. *
  992. *                    Generate left vectors bidiagonalizing R in WORK(IR)
  993. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  994. *
  995.                      CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  996.      $                            WORK( ITAUQ ), WORK( IWORK ),
  997.      $                            LWORK-IWORK+1, IERR )
  998.                      IWORK = IE + N
  999. *
  1000. *                    Perform bidiagonal QR iteration, computing left
  1001. *                    singular vectors of R in WORK(IR)
  1002. *                    (Workspace: need N*N+BDSPAC)
  1003. *
  1004.                      CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
  1005.      $                            1, WORK( IR ), LDWRKR, DUM, 1,
  1006.      $                            WORK( IWORK ), INFO )
  1007. *
  1008. *                    Multiply Q in A by left singular vectors of R in
  1009. *                    WORK(IR), storing result in U
  1010. *                    (Workspace: need N*N)
  1011. *
  1012.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
  1013.      $                           WORK( IR ), LDWRKR, ZERO, U, LDU )
  1014. *
  1015.                   ELSE
  1016. *
  1017. *                    Insufficient workspace for a fast algorithm
  1018. *
  1019.                      ITAU = 1
  1020.                      IWORK = ITAU + N
  1021. *
  1022. *                    Compute A=Q*R, copying result to U
  1023. *                    (Workspace: need 2*N, prefer N+N*NB)
  1024. *
  1025.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1026.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1027.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1028. *
  1029. *                    Generate Q in U
  1030. *                    (Workspace: need 2*N, prefer N+N*NB)
  1031. *
  1032.                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
  1033.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1034.                      IE = ITAU
  1035.                      ITAUQ = IE + N
  1036.                      ITAUP = ITAUQ + N
  1037.                      IWORK = ITAUP + N
  1038. *
  1039. *                    Zero out below R in A
  1040. *
  1041.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
  1042.      $                            LDA )
  1043. *
  1044. *                    Bidiagonalize R in A
  1045. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1046. *
  1047.                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
  1048.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1049.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1050. *
  1051. *                    Multiply Q in U by left vectors bidiagonalizing R
  1052. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1053. *
  1054.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1055.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1056.      $                            LWORK-IWORK+1, IERR )
  1057.                      IWORK = IE + N
  1058. *
  1059. *                    Perform bidiagonal QR iteration, computing left
  1060. *                    singular vectors of A in U
  1061. *                    (Workspace: need BDSPAC)
  1062. *
  1063.                      CALL DBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM,
  1064.      $                            1, U, LDU, DUM, 1, WORK( IWORK ),
  1065.      $                            INFO )
  1066. *
  1067.                   END IF
  1068. *
  1069.                ELSE IF( WNTVO ) THEN
  1070. *
  1071. *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
  1072. *                 N left singular vectors to be computed in U and
  1073. *                 N right singular vectors to be overwritten on A
  1074. *
  1075.                   IF( LWORK.GE.2*N*N+MAX( 4*N, BDSPAC ) ) THEN
  1076. *
  1077. *                    Sufficient workspace for a fast algorithm
  1078. *
  1079.                      IU = 1
  1080.                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
  1081. *
  1082. *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
  1083. *
  1084.                         LDWRKU = LDA
  1085.                         IR = IU + LDWRKU*N
  1086.                         LDWRKR = LDA
  1087.                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
  1088. *
  1089. *                       WORK(IU) is LDA by N and WORK(IR) is N by N
  1090. *
  1091.                         LDWRKU = LDA
  1092.                         IR = IU + LDWRKU*N
  1093.                         LDWRKR = N
  1094.                      ELSE
  1095. *
  1096. *                       WORK(IU) is N by N and WORK(IR) is N by N
  1097. *
  1098.                         LDWRKU = N
  1099.                         IR = IU + LDWRKU*N
  1100.                         LDWRKR = N
  1101.                      END IF
  1102.                      ITAU = IR + LDWRKR*N
  1103.                      IWORK = ITAU + N
  1104. *
  1105. *                    Compute A=Q*R
  1106. *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1107. *
  1108.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1109.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1110. *
  1111. *                    Copy R to WORK(IU), zeroing out below it
  1112. *
  1113.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1114.      $                            LDWRKU )
  1115.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  1116.      $                            WORK( IU+1 ), LDWRKU )
  1117. *
  1118. *                    Generate Q in A
  1119. *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1120. *
  1121.                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  1122.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1123.                      IE = ITAU
  1124.                      ITAUQ = IE + N
  1125.                      ITAUP = ITAUQ + N
  1126.                      IWORK = ITAUP + N
  1127. *
  1128. *                    Bidiagonalize R in WORK(IU), copying result to
  1129. *                    WORK(IR)
  1130. *                    (Workspace: need 2*N*N+4*N,
  1131. *                                prefer 2*N*N+3*N+2*N*NB)
  1132. *
  1133.                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1134.      $                            WORK( IE ), WORK( ITAUQ ),
  1135.      $                            WORK( ITAUP ), WORK( IWORK ),
  1136.      $                            LWORK-IWORK+1, IERR )
  1137.                      CALL DLACPY( 'U', N, N, WORK( IU ), LDWRKU,
  1138.      $                            WORK( IR ), LDWRKR )
  1139. *
  1140. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1141. *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
  1142. *
  1143.                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1144.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1145.      $                            LWORK-IWORK+1, IERR )
  1146. *
  1147. *                    Generate right bidiagonalizing vectors in WORK(IR)
  1148. *                    (Workspace: need 2*N*N+4*N-1,
  1149. *                                prefer 2*N*N+3*N+(N-1)*NB)
  1150. *
  1151.                      CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
  1152.      $                            WORK( ITAUP ), WORK( IWORK ),
  1153.      $                            LWORK-IWORK+1, IERR )
  1154.                      IWORK = IE + N
  1155. *
  1156. *                    Perform bidiagonal QR iteration, computing left
  1157. *                    singular vectors of R in WORK(IU) and computing
  1158. *                    right singular vectors of R in WORK(IR)
  1159. *                    (Workspace: need 2*N*N+BDSPAC)
  1160. *
  1161.                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
  1162.      $                            WORK( IR ), LDWRKR, WORK( IU ),
  1163.      $                            LDWRKU, DUM, 1, WORK( IWORK ), INFO )
  1164. *
  1165. *                    Multiply Q in A by left singular vectors of R in
  1166. *                    WORK(IU), storing result in U
  1167. *                    (Workspace: need N*N)
  1168. *
  1169.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
  1170.      $                           WORK( IU ), LDWRKU, ZERO, U, LDU )
  1171. *
  1172. *                    Copy right singular vectors of R to A
  1173. *                    (Workspace: need N*N)
  1174. *
  1175.                      CALL DLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
  1176.      $                            LDA )
  1177. *
  1178.                   ELSE
  1179. *
  1180. *                    Insufficient workspace for a fast algorithm
  1181. *
  1182.                      ITAU = 1
  1183.                      IWORK = ITAU + N
  1184. *
  1185. *                    Compute A=Q*R, copying result to U
  1186. *                    (Workspace: need 2*N, prefer N+N*NB)
  1187. *
  1188.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1189.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1190.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1191. *
  1192. *                    Generate Q in U
  1193. *                    (Workspace: need 2*N, prefer N+N*NB)
  1194. *
  1195.                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
  1196.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1197.                      IE = ITAU
  1198.                      ITAUQ = IE + N
  1199.                      ITAUP = ITAUQ + N
  1200.                      IWORK = ITAUP + N
  1201. *
  1202. *                    Zero out below R in A
  1203. *
  1204.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
  1205.      $                            LDA )
  1206. *
  1207. *                    Bidiagonalize R in A
  1208. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1209. *
  1210.                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
  1211.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1212.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1213. *
  1214. *                    Multiply Q in U by left vectors bidiagonalizing R
  1215. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1216. *
  1217.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1218.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1219.      $                            LWORK-IWORK+1, IERR )
  1220. *
  1221. *                    Generate right vectors bidiagonalizing R in A
  1222. *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1223. *
  1224.                      CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  1225.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1226.                      IWORK = IE + N
  1227. *
  1228. *                    Perform bidiagonal QR iteration, computing left
  1229. *                    singular vectors of A in U and computing right
  1230. *                    singular vectors of A in A
  1231. *                    (Workspace: need BDSPAC)
  1232. *
  1233.                      CALL DBDSQR( 'U', N, N, M, 0, S, WORK( IE ), A,
  1234.      $                            LDA, U, LDU, DUM, 1, WORK( IWORK ),
  1235.      $                            INFO )
  1236. *
  1237.                   END IF
  1238. *
  1239.                ELSE IF( WNTVAS ) THEN
  1240. *
  1241. *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
  1242. *                         or 'A')
  1243. *                 N left singular vectors to be computed in U and
  1244. *                 N right singular vectors to be computed in VT
  1245. *
  1246.                   IF( LWORK.GE.N*N+MAX( 4*N, BDSPAC ) ) THEN
  1247. *
  1248. *                    Sufficient workspace for a fast algorithm
  1249. *
  1250.                      IU = 1
  1251.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1252. *
  1253. *                       WORK(IU) is LDA by N
  1254. *
  1255.                         LDWRKU = LDA
  1256.                      ELSE
  1257. *
  1258. *                       WORK(IU) is N by N
  1259. *
  1260.                         LDWRKU = N
  1261.                      END IF
  1262.                      ITAU = IU + LDWRKU*N
  1263.                      IWORK = ITAU + N
  1264. *
  1265. *                    Compute A=Q*R
  1266. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  1267. *
  1268.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1269.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1270. *
  1271. *                    Copy R to WORK(IU), zeroing out below it
  1272. *
  1273.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1274.      $                            LDWRKU )
  1275.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  1276.      $                            WORK( IU+1 ), LDWRKU )
  1277. *
  1278. *                    Generate Q in A
  1279. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  1280. *
  1281.                      CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
  1282.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1283.                      IE = ITAU
  1284.                      ITAUQ = IE + N
  1285.                      ITAUP = ITAUQ + N
  1286.                      IWORK = ITAUP + N
  1287. *
  1288. *                    Bidiagonalize R in WORK(IU), copying result to VT
  1289. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  1290. *
  1291.                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1292.      $                            WORK( IE ), WORK( ITAUQ ),
  1293.      $                            WORK( ITAUP ), WORK( IWORK ),
  1294.      $                            LWORK-IWORK+1, IERR )
  1295.                      CALL DLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
  1296.      $                            LDVT )
  1297. *
  1298. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1299. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  1300. *
  1301.                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1302.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1303.      $                            LWORK-IWORK+1, IERR )
  1304. *
  1305. *                    Generate right bidiagonalizing vectors in VT
  1306. *                    (Workspace: need N*N+4*N-1,
  1307. *                                prefer N*N+3*N+(N-1)*NB)
  1308. *
  1309.                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1310.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1311.                      IWORK = IE + N
  1312. *
  1313. *                    Perform bidiagonal QR iteration, computing left
  1314. *                    singular vectors of R in WORK(IU) and computing
  1315. *                    right singular vectors of R in VT
  1316. *                    (Workspace: need N*N+BDSPAC)
  1317. *
  1318.                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
  1319.      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
  1320.      $                            WORK( IWORK ), INFO )
  1321. *
  1322. *                    Multiply Q in A by left singular vectors of R in
  1323. *                    WORK(IU), storing result in U
  1324. *                    (Workspace: need N*N)
  1325. *
  1326.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA,
  1327.      $                           WORK( IU ), LDWRKU, ZERO, U, LDU )
  1328. *
  1329.                   ELSE
  1330. *
  1331. *                    Insufficient workspace for a fast algorithm
  1332. *
  1333.                      ITAU = 1
  1334.                      IWORK = ITAU + N
  1335. *
  1336. *                    Compute A=Q*R, copying result to U
  1337. *                    (Workspace: need 2*N, prefer N+N*NB)
  1338. *
  1339.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1340.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1341.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1342. *
  1343. *                    Generate Q in U
  1344. *                    (Workspace: need 2*N, prefer N+N*NB)
  1345. *
  1346.                      CALL DORGQR( M, N, N, U, LDU, WORK( ITAU ),
  1347.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1348. *
  1349. *                    Copy R to VT, zeroing out below it
  1350. *
  1351.                      CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
  1352.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, VT( 2, 1 ),
  1353.      $                            LDVT )
  1354.                      IE = ITAU
  1355.                      ITAUQ = IE + N
  1356.                      ITAUP = ITAUQ + N
  1357.                      IWORK = ITAUP + N
  1358. *
  1359. *                    Bidiagonalize R in VT
  1360. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1361. *
  1362.                      CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
  1363.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1364.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1365. *
  1366. *                    Multiply Q in U by left bidiagonalizing vectors
  1367. *                    in VT
  1368. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1369. *
  1370.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  1371.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1372.      $                            LWORK-IWORK+1, IERR )
  1373. *
  1374. *                    Generate right bidiagonalizing vectors in VT
  1375. *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1376. *
  1377.                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1378.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1379.                      IWORK = IE + N
  1380. *
  1381. *                    Perform bidiagonal QR iteration, computing left
  1382. *                    singular vectors of A in U and computing right
  1383. *                    singular vectors of A in VT
  1384. *                    (Workspace: need BDSPAC)
  1385. *
  1386.                      CALL DBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT,
  1387.      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
  1388.      $                            INFO )
  1389. *
  1390.                   END IF
  1391. *
  1392.                END IF
  1393. *
  1394.             ELSE IF( WNTUA ) THEN
  1395. *
  1396.                IF( WNTVN ) THEN
  1397. *
  1398. *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
  1399. *                 M left singular vectors to be computed in U and
  1400. *                 no right singular vectors to be computed
  1401. *
  1402.                   IF( LWORK.GE.N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
  1403. *
  1404. *                    Sufficient workspace for a fast algorithm
  1405. *
  1406.                      IR = 1
  1407.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1408. *
  1409. *                       WORK(IR) is LDA by N
  1410. *
  1411.                         LDWRKR = LDA
  1412.                      ELSE
  1413. *
  1414. *                       WORK(IR) is N by N
  1415. *
  1416.                         LDWRKR = N
  1417.                      END IF
  1418.                      ITAU = IR + LDWRKR*N
  1419.                      IWORK = ITAU + N
  1420. *
  1421. *                    Compute A=Q*R, copying result to U
  1422. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  1423. *
  1424.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1425.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1426.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1427. *
  1428. *                    Copy R to WORK(IR), zeroing out below it
  1429. *
  1430.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ),
  1431.      $                            LDWRKR )
  1432.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  1433.      $                            WORK( IR+1 ), LDWRKR )
  1434. *
  1435. *                    Generate Q in U
  1436. *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
  1437. *
  1438.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1439.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1440.                      IE = ITAU
  1441.                      ITAUQ = IE + N
  1442.                      ITAUP = ITAUQ + N
  1443.                      IWORK = ITAUP + N
  1444. *
  1445. *                    Bidiagonalize R in WORK(IR)
  1446. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  1447. *
  1448.                      CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S,
  1449.      $                            WORK( IE ), WORK( ITAUQ ),
  1450.      $                            WORK( ITAUP ), WORK( IWORK ),
  1451.      $                            LWORK-IWORK+1, IERR )
  1452. *
  1453. *                    Generate left bidiagonalizing vectors in WORK(IR)
  1454. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  1455. *
  1456.                      CALL DORGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
  1457.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1458.      $                            LWORK-IWORK+1, IERR )
  1459.                      IWORK = IE + N
  1460. *
  1461. *                    Perform bidiagonal QR iteration, computing left
  1462. *                    singular vectors of R in WORK(IR)
  1463. *                    (Workspace: need N*N+BDSPAC)
  1464. *
  1465.                      CALL DBDSQR( 'U', N, 0, N, 0, S, WORK( IE ), DUM,
  1466.      $                            1, WORK( IR ), LDWRKR, DUM, 1,
  1467.      $                            WORK( IWORK ), INFO )
  1468. *
  1469. *                    Multiply Q in U by left singular vectors of R in
  1470. *                    WORK(IR), storing result in A
  1471. *                    (Workspace: need N*N)
  1472. *
  1473.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
  1474.      $                           WORK( IR ), LDWRKR, ZERO, A, LDA )
  1475. *
  1476. *                    Copy left singular vectors of A from A to U
  1477. *
  1478.                      CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
  1479. *
  1480.                   ELSE
  1481. *
  1482. *                    Insufficient workspace for a fast algorithm
  1483. *
  1484.                      ITAU = 1
  1485.                      IWORK = ITAU + N
  1486. *
  1487. *                    Compute A=Q*R, copying result to U
  1488. *                    (Workspace: need 2*N, prefer N+N*NB)
  1489. *
  1490.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1491.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1492.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1493. *
  1494. *                    Generate Q in U
  1495. *                    (Workspace: need N+M, prefer N+M*NB)
  1496. *
  1497.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1498.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1499.                      IE = ITAU
  1500.                      ITAUQ = IE + N
  1501.                      ITAUP = ITAUQ + N
  1502.                      IWORK = ITAUP + N
  1503. *
  1504. *                    Zero out below R in A
  1505. *
  1506.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
  1507.      $                            LDA )
  1508. *
  1509. *                    Bidiagonalize R in A
  1510. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1511. *
  1512.                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
  1513.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1514.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1515. *
  1516. *                    Multiply Q in U by left bidiagonalizing vectors
  1517. *                    in A
  1518. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1519. *
  1520.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1521.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1522.      $                            LWORK-IWORK+1, IERR )
  1523.                      IWORK = IE + N
  1524. *
  1525. *                    Perform bidiagonal QR iteration, computing left
  1526. *                    singular vectors of A in U
  1527. *                    (Workspace: need BDSPAC)
  1528. *
  1529.                      CALL DBDSQR( 'U', N, 0, M, 0, S, WORK( IE ), DUM,
  1530.      $                            1, U, LDU, DUM, 1, WORK( IWORK ),
  1531.      $                            INFO )
  1532. *
  1533.                   END IF
  1534. *
  1535.                ELSE IF( WNTVO ) THEN
  1536. *
  1537. *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
  1538. *                 M left singular vectors to be computed in U and
  1539. *                 N right singular vectors to be overwritten on A
  1540. *
  1541.                   IF( LWORK.GE.2*N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
  1542. *
  1543. *                    Sufficient workspace for a fast algorithm
  1544. *
  1545.                      IU = 1
  1546.                      IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
  1547. *
  1548. *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
  1549. *
  1550.                         LDWRKU = LDA
  1551.                         IR = IU + LDWRKU*N
  1552.                         LDWRKR = LDA
  1553.                      ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
  1554. *
  1555. *                       WORK(IU) is LDA by N and WORK(IR) is N by N
  1556. *
  1557.                         LDWRKU = LDA
  1558.                         IR = IU + LDWRKU*N
  1559.                         LDWRKR = N
  1560.                      ELSE
  1561. *
  1562. *                       WORK(IU) is N by N and WORK(IR) is N by N
  1563. *
  1564.                         LDWRKU = N
  1565.                         IR = IU + LDWRKU*N
  1566.                         LDWRKR = N
  1567.                      END IF
  1568.                      ITAU = IR + LDWRKR*N
  1569.                      IWORK = ITAU + N
  1570. *
  1571. *                    Compute A=Q*R, copying result to U
  1572. *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
  1573. *
  1574.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1575.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1576.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1577. *
  1578. *                    Generate Q in U
  1579. *                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
  1580. *
  1581.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1582.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1583. *
  1584. *                    Copy R to WORK(IU), zeroing out below it
  1585. *
  1586.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1587.      $                            LDWRKU )
  1588.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  1589.      $                            WORK( IU+1 ), LDWRKU )
  1590.                      IE = ITAU
  1591.                      ITAUQ = IE + N
  1592.                      ITAUP = ITAUQ + N
  1593.                      IWORK = ITAUP + N
  1594. *
  1595. *                    Bidiagonalize R in WORK(IU), copying result to
  1596. *                    WORK(IR)
  1597. *                    (Workspace: need 2*N*N+4*N,
  1598. *                                prefer 2*N*N+3*N+2*N*NB)
  1599. *
  1600.                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1601.      $                            WORK( IE ), WORK( ITAUQ ),
  1602.      $                            WORK( ITAUP ), WORK( IWORK ),
  1603.      $                            LWORK-IWORK+1, IERR )
  1604.                      CALL DLACPY( 'U', N, N, WORK( IU ), LDWRKU,
  1605.      $                            WORK( IR ), LDWRKR )
  1606. *
  1607. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1608. *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
  1609. *
  1610.                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1611.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1612.      $                            LWORK-IWORK+1, IERR )
  1613. *
  1614. *                    Generate right bidiagonalizing vectors in WORK(IR)
  1615. *                    (Workspace: need 2*N*N+4*N-1,
  1616. *                                prefer 2*N*N+3*N+(N-1)*NB)
  1617. *
  1618.                      CALL DORGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
  1619.      $                            WORK( ITAUP ), WORK( IWORK ),
  1620.      $                            LWORK-IWORK+1, IERR )
  1621.                      IWORK = IE + N
  1622. *
  1623. *                    Perform bidiagonal QR iteration, computing left
  1624. *                    singular vectors of R in WORK(IU) and computing
  1625. *                    right singular vectors of R in WORK(IR)
  1626. *                    (Workspace: need 2*N*N+BDSPAC)
  1627. *
  1628.                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ),
  1629.      $                            WORK( IR ), LDWRKR, WORK( IU ),
  1630.      $                            LDWRKU, DUM, 1, WORK( IWORK ), INFO )
  1631. *
  1632. *                    Multiply Q in U by left singular vectors of R in
  1633. *                    WORK(IU), storing result in A
  1634. *                    (Workspace: need N*N)
  1635. *
  1636.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
  1637.      $                           WORK( IU ), LDWRKU, ZERO, A, LDA )
  1638. *
  1639. *                    Copy left singular vectors of A from A to U
  1640. *
  1641.                      CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
  1642. *
  1643. *                    Copy right singular vectors of R from WORK(IR) to A
  1644. *
  1645.                      CALL DLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
  1646.      $                            LDA )
  1647. *
  1648.                   ELSE
  1649. *
  1650. *                    Insufficient workspace for a fast algorithm
  1651. *
  1652.                      ITAU = 1
  1653.                      IWORK = ITAU + N
  1654. *
  1655. *                    Compute A=Q*R, copying result to U
  1656. *                    (Workspace: need 2*N, prefer N+N*NB)
  1657. *
  1658.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1659.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1660.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1661. *
  1662. *                    Generate Q in U
  1663. *                    (Workspace: need N+M, prefer N+M*NB)
  1664. *
  1665.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1666.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1667.                      IE = ITAU
  1668.                      ITAUQ = IE + N
  1669.                      ITAUP = ITAUQ + N
  1670.                      IWORK = ITAUP + N
  1671. *
  1672. *                    Zero out below R in A
  1673. *
  1674.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ),
  1675.      $                            LDA )
  1676. *
  1677. *                    Bidiagonalize R in A
  1678. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1679. *
  1680.                      CALL DGEBRD( N, N, A, LDA, S, WORK( IE ),
  1681.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1682.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1683. *
  1684. *                    Multiply Q in U by left bidiagonalizing vectors
  1685. *                    in A
  1686. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1687. *
  1688.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
  1689.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1690.      $                            LWORK-IWORK+1, IERR )
  1691. *
  1692. *                    Generate right bidiagonalizing vectors in A
  1693. *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1694. *
  1695.                      CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  1696.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1697.                      IWORK = IE + N
  1698. *
  1699. *                    Perform bidiagonal QR iteration, computing left
  1700. *                    singular vectors of A in U and computing right
  1701. *                    singular vectors of A in A
  1702. *                    (Workspace: need BDSPAC)
  1703. *
  1704.                      CALL DBDSQR( 'U', N, N, M, 0, S, WORK( IE ), A,
  1705.      $                            LDA, U, LDU, DUM, 1, WORK( IWORK ),
  1706.      $                            INFO )
  1707. *
  1708.                   END IF
  1709. *
  1710.                ELSE IF( WNTVAS ) THEN
  1711. *
  1712. *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
  1713. *                         or 'A')
  1714. *                 M left singular vectors to be computed in U and
  1715. *                 N right singular vectors to be computed in VT
  1716. *
  1717.                   IF( LWORK.GE.N*N+MAX( N+M, 4*N, BDSPAC ) ) THEN
  1718. *
  1719. *                    Sufficient workspace for a fast algorithm
  1720. *
  1721.                      IU = 1
  1722.                      IF( LWORK.GE.WRKBL+LDA*N ) THEN
  1723. *
  1724. *                       WORK(IU) is LDA by N
  1725. *
  1726.                         LDWRKU = LDA
  1727.                      ELSE
  1728. *
  1729. *                       WORK(IU) is N by N
  1730. *
  1731.                         LDWRKU = N
  1732.                      END IF
  1733.                      ITAU = IU + LDWRKU*N
  1734.                      IWORK = ITAU + N
  1735. *
  1736. *                    Compute A=Q*R, copying result to U
  1737. *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
  1738. *
  1739.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1740.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1741.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1742. *
  1743. *                    Generate Q in U
  1744. *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
  1745. *
  1746.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1747.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1748. *
  1749. *                    Copy R to WORK(IU), zeroing out below it
  1750. *
  1751.                      CALL DLACPY( 'U', N, N, A, LDA, WORK( IU ),
  1752.      $                            LDWRKU )
  1753.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO,
  1754.      $                            WORK( IU+1 ), LDWRKU )
  1755.                      IE = ITAU
  1756.                      ITAUQ = IE + N
  1757.                      ITAUP = ITAUQ + N
  1758.                      IWORK = ITAUP + N
  1759. *
  1760. *                    Bidiagonalize R in WORK(IU), copying result to VT
  1761. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
  1762. *
  1763.                      CALL DGEBRD( N, N, WORK( IU ), LDWRKU, S,
  1764.      $                            WORK( IE ), WORK( ITAUQ ),
  1765.      $                            WORK( ITAUP ), WORK( IWORK ),
  1766.      $                            LWORK-IWORK+1, IERR )
  1767.                      CALL DLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
  1768.      $                            LDVT )
  1769. *
  1770. *                    Generate left bidiagonalizing vectors in WORK(IU)
  1771. *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
  1772. *
  1773.                      CALL DORGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
  1774.      $                            WORK( ITAUQ ), WORK( IWORK ),
  1775.      $                            LWORK-IWORK+1, IERR )
  1776. *
  1777. *                    Generate right bidiagonalizing vectors in VT
  1778. *                    (Workspace: need N*N+4*N-1,
  1779. *                                prefer N*N+3*N+(N-1)*NB)
  1780. *
  1781.                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1782.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1783.                      IWORK = IE + N
  1784. *
  1785. *                    Perform bidiagonal QR iteration, computing left
  1786. *                    singular vectors of R in WORK(IU) and computing
  1787. *                    right singular vectors of R in VT
  1788. *                    (Workspace: need N*N+BDSPAC)
  1789. *
  1790.                      CALL DBDSQR( 'U', N, N, N, 0, S, WORK( IE ), VT,
  1791.      $                            LDVT, WORK( IU ), LDWRKU, DUM, 1,
  1792.      $                            WORK( IWORK ), INFO )
  1793. *
  1794. *                    Multiply Q in U by left singular vectors of R in
  1795. *                    WORK(IU), storing result in A
  1796. *                    (Workspace: need N*N)
  1797. *
  1798.                      CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU,
  1799.      $                           WORK( IU ), LDWRKU, ZERO, A, LDA )
  1800. *
  1801. *                    Copy left singular vectors of A from A to U
  1802. *
  1803.                      CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
  1804. *
  1805.                   ELSE
  1806. *
  1807. *                    Insufficient workspace for a fast algorithm
  1808. *
  1809.                      ITAU = 1
  1810.                      IWORK = ITAU + N
  1811. *
  1812. *                    Compute A=Q*R, copying result to U
  1813. *                    (Workspace: need 2*N, prefer N+N*NB)
  1814. *
  1815.                      CALL DGEQRF( M, N, A, LDA, WORK( ITAU ),
  1816.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1817.                      CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1818. *
  1819. *                    Generate Q in U
  1820. *                    (Workspace: need N+M, prefer N+M*NB)
  1821. *
  1822.                      CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
  1823.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1824. *
  1825. *                    Copy R from A to VT, zeroing out below it
  1826. *
  1827.                      CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
  1828.                      CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, VT( 2, 1 ),
  1829.      $                            LDVT )
  1830.                      IE = ITAU
  1831.                      ITAUQ = IE + N
  1832.                      ITAUP = ITAUQ + N
  1833.                      IWORK = ITAUP + N
  1834. *
  1835. *                    Bidiagonalize R in VT
  1836. *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
  1837. *
  1838.                      CALL DGEBRD( N, N, VT, LDVT, S, WORK( IE ),
  1839.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  1840.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1841. *
  1842. *                    Multiply Q in U by left bidiagonalizing vectors
  1843. *                    in VT
  1844. *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
  1845. *
  1846.                      CALL DORMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
  1847.      $                            WORK( ITAUQ ), U, LDU, WORK( IWORK ),
  1848.      $                            LWORK-IWORK+1, IERR )
  1849. *
  1850. *                    Generate right bidiagonalizing vectors in VT
  1851. *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1852. *
  1853.                      CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1854.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  1855.                      IWORK = IE + N
  1856. *
  1857. *                    Perform bidiagonal QR iteration, computing left
  1858. *                    singular vectors of A in U and computing right
  1859. *                    singular vectors of A in VT
  1860. *                    (Workspace: need BDSPAC)
  1861. *
  1862.                      CALL DBDSQR( 'U', N, N, M, 0, S, WORK( IE ), VT,
  1863.      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
  1864.      $                            INFO )
  1865. *
  1866.                   END IF
  1867. *
  1868.                END IF
  1869. *
  1870.             END IF
  1871. *
  1872.          ELSE
  1873. *
  1874. *           M .LT. MNTHR
  1875. *
  1876. *           Path 10 (M at least N, but not much larger)
  1877. *           Reduce to bidiagonal form without QR decomposition
  1878. *
  1879.             IE = 1
  1880.             ITAUQ = IE + N
  1881.             ITAUP = ITAUQ + N
  1882.             IWORK = ITAUP + N
  1883. *
  1884. *           Bidiagonalize A
  1885. *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
  1886. *
  1887.             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  1888.      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  1889.      $                   IERR )
  1890.             IF( WNTUAS ) THEN
  1891. *
  1892. *              If left singular vectors desired in U, copy result to U
  1893. *              and generate left bidiagonalizing vectors in U
  1894. *              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
  1895. *
  1896.                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
  1897.                IF( WNTUS )
  1898.      $            NCU = N
  1899.                IF( WNTUA )
  1900.      $            NCU = M
  1901.                CALL DORGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
  1902.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  1903.             END IF
  1904.             IF( WNTVAS ) THEN
  1905. *
  1906. *              If right singular vectors desired in VT, copy result to
  1907. *              VT and generate right bidiagonalizing vectors in VT
  1908. *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1909. *
  1910.                CALL DLACPY( 'U', N, N, A, LDA, VT, LDVT )
  1911.                CALL DORGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
  1912.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  1913.             END IF
  1914.             IF( WNTUO ) THEN
  1915. *
  1916. *              If left singular vectors desired in A, generate left
  1917. *              bidiagonalizing vectors in A
  1918. *              (Workspace: need 4*N, prefer 3*N+N*NB)
  1919. *
  1920.                CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
  1921.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  1922.             END IF
  1923.             IF( WNTVO ) THEN
  1924. *
  1925. *              If right singular vectors desired in A, generate right
  1926. *              bidiagonalizing vectors in A
  1927. *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
  1928. *
  1929.                CALL DORGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
  1930.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  1931.             END IF
  1932.             IWORK = IE + N
  1933.             IF( WNTUAS .OR. WNTUO )
  1934.      $         NRU = M
  1935.             IF( WNTUN )
  1936.      $         NRU = 0
  1937.             IF( WNTVAS .OR. WNTVO )
  1938.      $         NCVT = N
  1939.             IF( WNTVN )
  1940.      $         NCVT = 0
  1941.             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
  1942. *
  1943. *              Perform bidiagonal QR iteration, if desired, computing
  1944. *              left singular vectors in U and computing right singular
  1945. *              vectors in VT
  1946. *              (Workspace: need BDSPAC)
  1947. *
  1948.                CALL DBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), VT,
  1949.      $                      LDVT, U, LDU, DUM, 1, WORK( IWORK ), INFO )
  1950.             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
  1951. *
  1952. *              Perform bidiagonal QR iteration, if desired, computing
  1953. *              left singular vectors in U and computing right singular
  1954. *              vectors in A
  1955. *              (Workspace: need BDSPAC)
  1956. *
  1957.                CALL DBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), A, LDA,
  1958.      $                      U, LDU, DUM, 1, WORK( IWORK ), INFO )
  1959.             ELSE
  1960. *
  1961. *              Perform bidiagonal QR iteration, if desired, computing
  1962. *              left singular vectors in A and computing right singular
  1963. *              vectors in VT
  1964. *              (Workspace: need BDSPAC)
  1965. *
  1966.                CALL DBDSQR( 'U', N, NCVT, NRU, 0, S, WORK( IE ), VT,
  1967.      $                      LDVT, A, LDA, DUM, 1, WORK( IWORK ), INFO )
  1968.             END IF
  1969. *
  1970.          END IF
  1971. *
  1972.       ELSE
  1973. *
  1974. *        A has more columns than rows. If A has sufficiently more
  1975. *        columns than rows, first reduce using the LQ decomposition (if
  1976. *        sufficient workspace available)
  1977. *
  1978.          IF( N.GE.MNTHR ) THEN
  1979. *
  1980.             IF( WNTVN ) THEN
  1981. *
  1982. *              Path 1t(N much larger than M, JOBVT='N')
  1983. *              No right singular vectors to be computed
  1984. *
  1985.                ITAU = 1
  1986.                IWORK = ITAU + M
  1987. *
  1988. *              Compute A=L*Q
  1989. *              (Workspace: need 2*M, prefer M+M*NB)
  1990. *
  1991.                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
  1992.      $                      LWORK-IWORK+1, IERR )
  1993. *
  1994. *              Zero out above L
  1995. *
  1996.                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
  1997.                IE = 1
  1998.                ITAUQ = IE + M
  1999.                ITAUP = ITAUQ + M
  2000.                IWORK = ITAUP + M
  2001. *
  2002. *              Bidiagonalize L in A
  2003. *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2004. *
  2005.                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  2006.      $                      WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  2007.      $                      IERR )
  2008.                IF( WNTUO .OR. WNTUAS ) THEN
  2009. *
  2010. *                 If left singular vectors desired, generate Q
  2011. *                 (Workspace: need 4*M, prefer 3*M+M*NB)
  2012. *
  2013.                   CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  2014.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2015.                END IF
  2016.                IWORK = IE + M
  2017.                NRU = 0
  2018.                IF( WNTUO .OR. WNTUAS )
  2019.      $            NRU = M
  2020. *
  2021. *              Perform bidiagonal QR iteration, computing left singular
  2022. *              vectors of A in A if desired
  2023. *              (Workspace: need BDSPAC)
  2024. *
  2025.                CALL DBDSQR( 'U', M, 0, NRU, 0, S, WORK( IE ), DUM, 1, A,
  2026.      $                      LDA, DUM, 1, WORK( IWORK ), INFO )
  2027. *
  2028. *              If left singular vectors desired in U, copy them there
  2029. *
  2030.                IF( WNTUAS )
  2031.      $            CALL DLACPY( 'F', M, M, A, LDA, U, LDU )
  2032. *
  2033.             ELSE IF( WNTVO .AND. WNTUN ) THEN
  2034. *
  2035. *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
  2036. *              M right singular vectors to be overwritten on A and
  2037. *              no left singular vectors to be computed
  2038. *
  2039.                IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
  2040. *
  2041. *                 Sufficient workspace for a fast algorithm
  2042. *
  2043.                   IR = 1
  2044.                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
  2045. *
  2046. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
  2047. *
  2048.                      LDWRKU = LDA
  2049.                      CHUNK = N
  2050.                      LDWRKR = LDA
  2051.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
  2052. *
  2053. *                    WORK(IU) is LDA by N and WORK(IR) is M by M
  2054. *
  2055.                      LDWRKU = LDA
  2056.                      CHUNK = N
  2057.                      LDWRKR = M
  2058.                   ELSE
  2059. *
  2060. *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
  2061. *
  2062.                      LDWRKU = M
  2063.                      CHUNK = ( LWORK-M*M-M ) / M
  2064.                      LDWRKR = M
  2065.                   END IF
  2066.                   ITAU = IR + LDWRKR*M
  2067.                   IWORK = ITAU + M
  2068. *
  2069. *                 Compute A=L*Q
  2070. *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2071. *
  2072.                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2073.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2074. *
  2075. *                 Copy L to WORK(IR) and zero out above it
  2076. *
  2077.                   CALL DLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
  2078.                   CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2079.      $                         WORK( IR+LDWRKR ), LDWRKR )
  2080. *
  2081. *                 Generate Q in A
  2082. *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2083. *
  2084.                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2085.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2086.                   IE = ITAU
  2087.                   ITAUQ = IE + M
  2088.                   ITAUP = ITAUQ + M
  2089.                   IWORK = ITAUP + M
  2090. *
  2091. *                 Bidiagonalize L in WORK(IR)
  2092. *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  2093. *
  2094.                   CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S, WORK( IE ),
  2095.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2096.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2097. *
  2098. *                 Generate right vectors bidiagonalizing L
  2099. *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
  2100. *
  2101.                   CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2102.      $                         WORK( ITAUP ), WORK( IWORK ),
  2103.      $                         LWORK-IWORK+1, IERR )
  2104.                   IWORK = IE + M
  2105. *
  2106. *                 Perform bidiagonal QR iteration, computing right
  2107. *                 singular vectors of L in WORK(IR)
  2108. *                 (Workspace: need M*M+BDSPAC)
  2109. *
  2110.                   CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
  2111.      $                         WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
  2112.      $                         WORK( IWORK ), INFO )
  2113.                   IU = IE + M
  2114. *
  2115. *                 Multiply right singular vectors of L in WORK(IR) by Q
  2116. *                 in A, storing result in WORK(IU) and copying to A
  2117. *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)
  2118. *
  2119.                   DO 30 I = 1, N, CHUNK
  2120.                      BLK = MIN( N-I+1, CHUNK )
  2121.                      CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IR ),
  2122.      $                           LDWRKR, A( 1, I ), LDA, ZERO,
  2123.      $                           WORK( IU ), LDWRKU )
  2124.                      CALL DLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
  2125.      $                            A( 1, I ), LDA )
  2126.    30             CONTINUE
  2127. *
  2128.                ELSE
  2129. *
  2130. *                 Insufficient workspace for a fast algorithm
  2131. *
  2132.                   IE = 1
  2133.                   ITAUQ = IE + M
  2134.                   ITAUP = ITAUQ + M
  2135.                   IWORK = ITAUP + M
  2136. *
  2137. *                 Bidiagonalize A
  2138. *                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
  2139. *
  2140.                   CALL DGEBRD( M, N, A, LDA, S, WORK( IE ),
  2141.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2142.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2143. *
  2144. *                 Generate right vectors bidiagonalizing A
  2145. *                 (Workspace: need 4*M, prefer 3*M+M*NB)
  2146. *
  2147.                   CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
  2148.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2149.                   IWORK = IE + M
  2150. *
  2151. *                 Perform bidiagonal QR iteration, computing right
  2152. *                 singular vectors of A in A
  2153. *                 (Workspace: need BDSPAC)
  2154. *
  2155.                   CALL DBDSQR( 'L', M, N, 0, 0, S, WORK( IE ), A, LDA,
  2156.      $                         DUM, 1, DUM, 1, WORK( IWORK ), INFO )
  2157. *
  2158.                END IF
  2159. *
  2160.             ELSE IF( WNTVO .AND. WNTUAS ) THEN
  2161. *
  2162. *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
  2163. *              M right singular vectors to be overwritten on A and
  2164. *              M left singular vectors to be computed in U
  2165. *
  2166.                IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
  2167. *
  2168. *                 Sufficient workspace for a fast algorithm
  2169. *
  2170.                   IR = 1
  2171.                   IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+LDA*M ) THEN
  2172. *
  2173. *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
  2174. *
  2175.                      LDWRKU = LDA
  2176.                      CHUNK = N
  2177.                      LDWRKR = LDA
  2178.                   ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N+M )+M*M ) THEN
  2179. *
  2180. *                    WORK(IU) is LDA by N and WORK(IR) is M by M
  2181. *
  2182.                      LDWRKU = LDA
  2183.                      CHUNK = N
  2184.                      LDWRKR = M
  2185.                   ELSE
  2186. *
  2187. *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
  2188. *
  2189.                      LDWRKU = M
  2190.                      CHUNK = ( LWORK-M*M-M ) / M
  2191.                      LDWRKR = M
  2192.                   END IF
  2193.                   ITAU = IR + LDWRKR*M
  2194.                   IWORK = ITAU + M
  2195. *
  2196. *                 Compute A=L*Q
  2197. *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2198. *
  2199.                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2200.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2201. *
  2202. *                 Copy L to U, zeroing about above it
  2203. *
  2204.                   CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
  2205.                   CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
  2206.      $                         LDU )
  2207. *
  2208. *                 Generate Q in A
  2209. *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2210. *
  2211.                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2212.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2213.                   IE = ITAU
  2214.                   ITAUQ = IE + M
  2215.                   ITAUP = ITAUQ + M
  2216.                   IWORK = ITAUP + M
  2217. *
  2218. *                 Bidiagonalize L in U, copying result to WORK(IR)
  2219. *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  2220. *
  2221.                   CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
  2222.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2223.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2224.                   CALL DLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
  2225. *
  2226. *                 Generate right vectors bidiagonalizing L in WORK(IR)
  2227. *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
  2228. *
  2229.                   CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2230.      $                         WORK( ITAUP ), WORK( IWORK ),
  2231.      $                         LWORK-IWORK+1, IERR )
  2232. *
  2233. *                 Generate left vectors bidiagonalizing L in U
  2234. *                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
  2235. *
  2236.                   CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2237.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2238.                   IWORK = IE + M
  2239. *
  2240. *                 Perform bidiagonal QR iteration, computing left
  2241. *                 singular vectors of L in U, and computing right
  2242. *                 singular vectors of L in WORK(IR)
  2243. *                 (Workspace: need M*M+BDSPAC)
  2244. *
  2245.                   CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
  2246.      $                         WORK( IR ), LDWRKR, U, LDU, DUM, 1,
  2247.      $                         WORK( IWORK ), INFO )
  2248.                   IU = IE + M
  2249. *
  2250. *                 Multiply right singular vectors of L in WORK(IR) by Q
  2251. *                 in A, storing result in WORK(IU) and copying to A
  2252. *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M))
  2253. *
  2254.                   DO 40 I = 1, N, CHUNK
  2255.                      BLK = MIN( N-I+1, CHUNK )
  2256.                      CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IR ),
  2257.      $                           LDWRKR, A( 1, I ), LDA, ZERO,
  2258.      $                           WORK( IU ), LDWRKU )
  2259.                      CALL DLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
  2260.      $                            A( 1, I ), LDA )
  2261.    40             CONTINUE
  2262. *
  2263.                ELSE
  2264. *
  2265. *                 Insufficient workspace for a fast algorithm
  2266. *
  2267.                   ITAU = 1
  2268.                   IWORK = ITAU + M
  2269. *
  2270. *                 Compute A=L*Q
  2271. *                 (Workspace: need 2*M, prefer M+M*NB)
  2272. *
  2273.                   CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2274.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2275. *
  2276. *                 Copy L to U, zeroing out above it
  2277. *
  2278.                   CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
  2279.                   CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
  2280.      $                         LDU )
  2281. *
  2282. *                 Generate Q in A
  2283. *                 (Workspace: need 2*M, prefer M+M*NB)
  2284. *
  2285.                   CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2286.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2287.                   IE = ITAU
  2288.                   ITAUQ = IE + M
  2289.                   ITAUP = ITAUQ + M
  2290.                   IWORK = ITAUP + M
  2291. *
  2292. *                 Bidiagonalize L in U
  2293. *                 (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2294. *
  2295.                   CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
  2296.      $                         WORK( ITAUQ ), WORK( ITAUP ),
  2297.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2298. *
  2299. *                 Multiply right vectors bidiagonalizing L by Q in A
  2300. *                 (Workspace: need 3*M+N, prefer 3*M+N*NB)
  2301. *
  2302.                   CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
  2303.      $                         WORK( ITAUP ), A, LDA, WORK( IWORK ),
  2304.      $                         LWORK-IWORK+1, IERR )
  2305. *
  2306. *                 Generate left vectors bidiagonalizing L in U
  2307. *                 (Workspace: need 4*M, prefer 3*M+M*NB)
  2308. *
  2309.                   CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2310.      $                         WORK( IWORK ), LWORK-IWORK+1, IERR )
  2311.                   IWORK = IE + M
  2312. *
  2313. *                 Perform bidiagonal QR iteration, computing left
  2314. *                 singular vectors of A in U and computing right
  2315. *                 singular vectors of A in A
  2316. *                 (Workspace: need BDSPAC)
  2317. *
  2318.                   CALL DBDSQR( 'U', M, N, M, 0, S, WORK( IE ), A, LDA,
  2319.      $                         U, LDU, DUM, 1, WORK( IWORK ), INFO )
  2320. *
  2321.                END IF
  2322. *
  2323.             ELSE IF( WNTVS ) THEN
  2324. *
  2325.                IF( WNTUN ) THEN
  2326. *
  2327. *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
  2328. *                 M right singular vectors to be computed in VT and
  2329. *                 no left singular vectors to be computed
  2330. *
  2331.                   IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
  2332. *
  2333. *                    Sufficient workspace for a fast algorithm
  2334. *
  2335.                      IR = 1
  2336.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2337. *
  2338. *                       WORK(IR) is LDA by M
  2339. *
  2340.                         LDWRKR = LDA
  2341.                      ELSE
  2342. *
  2343. *                       WORK(IR) is M by M
  2344. *
  2345.                         LDWRKR = M
  2346.                      END IF
  2347.                      ITAU = IR + LDWRKR*M
  2348.                      IWORK = ITAU + M
  2349. *
  2350. *                    Compute A=L*Q
  2351. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2352. *
  2353.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2354.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2355. *
  2356. *                    Copy L to WORK(IR), zeroing out above it
  2357. *
  2358.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IR ),
  2359.      $                            LDWRKR )
  2360.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2361.      $                            WORK( IR+LDWRKR ), LDWRKR )
  2362. *
  2363. *                    Generate Q in A
  2364. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2365. *
  2366.                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2367.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2368.                      IE = ITAU
  2369.                      ITAUQ = IE + M
  2370.                      ITAUP = ITAUQ + M
  2371.                      IWORK = ITAUP + M
  2372. *
  2373. *                    Bidiagonalize L in WORK(IR)
  2374. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  2375. *
  2376.                      CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
  2377.      $                            WORK( IE ), WORK( ITAUQ ),
  2378.      $                            WORK( ITAUP ), WORK( IWORK ),
  2379.      $                            LWORK-IWORK+1, IERR )
  2380. *
  2381. *                    Generate right vectors bidiagonalizing L in
  2382. *                    WORK(IR)
  2383. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
  2384. *
  2385.                      CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2386.      $                            WORK( ITAUP ), WORK( IWORK ),
  2387.      $                            LWORK-IWORK+1, IERR )
  2388.                      IWORK = IE + M
  2389. *
  2390. *                    Perform bidiagonal QR iteration, computing right
  2391. *                    singular vectors of L in WORK(IR)
  2392. *                    (Workspace: need M*M+BDSPAC)
  2393. *
  2394.                      CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
  2395.      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
  2396.      $                            WORK( IWORK ), INFO )
  2397. *
  2398. *                    Multiply right singular vectors of L in WORK(IR) by
  2399. *                    Q in A, storing result in VT
  2400. *                    (Workspace: need M*M)
  2401. *
  2402.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IR ),
  2403.      $                           LDWRKR, A, LDA, ZERO, VT, LDVT )
  2404. *
  2405.                   ELSE
  2406. *
  2407. *                    Insufficient workspace for a fast algorithm
  2408. *
  2409.                      ITAU = 1
  2410.                      IWORK = ITAU + M
  2411. *
  2412. *                    Compute A=L*Q
  2413. *                    (Workspace: need 2*M, prefer M+M*NB)
  2414. *
  2415.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2416.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2417. *
  2418. *                    Copy result to VT
  2419. *
  2420.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2421. *
  2422. *                    Generate Q in VT
  2423. *                    (Workspace: need 2*M, prefer M+M*NB)
  2424. *
  2425.                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2426.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2427.                      IE = ITAU
  2428.                      ITAUQ = IE + M
  2429.                      ITAUP = ITAUQ + M
  2430.                      IWORK = ITAUP + M
  2431. *
  2432. *                    Zero out above L in A
  2433. *
  2434.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
  2435.      $                            LDA )
  2436. *
  2437. *                    Bidiagonalize L in A
  2438. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2439. *
  2440.                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
  2441.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2442.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2443. *
  2444. *                    Multiply right vectors bidiagonalizing L by Q in VT
  2445. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  2446. *
  2447.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
  2448.      $                            WORK( ITAUP ), VT, LDVT,
  2449.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2450.                      IWORK = IE + M
  2451. *
  2452. *                    Perform bidiagonal QR iteration, computing right
  2453. *                    singular vectors of A in VT
  2454. *                    (Workspace: need BDSPAC)
  2455. *
  2456.                      CALL DBDSQR( 'U', M, N, 0, 0, S, WORK( IE ), VT,
  2457.      $                            LDVT, DUM, 1, DUM, 1, WORK( IWORK ),
  2458.      $                            INFO )
  2459. *
  2460.                   END IF
  2461. *
  2462.                ELSE IF( WNTUO ) THEN
  2463. *
  2464. *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
  2465. *                 M right singular vectors to be computed in VT and
  2466. *                 M left singular vectors to be overwritten on A
  2467. *
  2468.                   IF( LWORK.GE.2*M*M+MAX( 4*M, BDSPAC ) ) THEN
  2469. *
  2470. *                    Sufficient workspace for a fast algorithm
  2471. *
  2472.                      IU = 1
  2473.                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
  2474. *
  2475. *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
  2476. *
  2477.                         LDWRKU = LDA
  2478.                         IR = IU + LDWRKU*M
  2479.                         LDWRKR = LDA
  2480.                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
  2481. *
  2482. *                       WORK(IU) is LDA by M and WORK(IR) is M by M
  2483. *
  2484.                         LDWRKU = LDA
  2485.                         IR = IU + LDWRKU*M
  2486.                         LDWRKR = M
  2487.                      ELSE
  2488. *
  2489. *                       WORK(IU) is M by M and WORK(IR) is M by M
  2490. *
  2491.                         LDWRKU = M
  2492.                         IR = IU + LDWRKU*M
  2493.                         LDWRKR = M
  2494.                      END IF
  2495.                      ITAU = IR + LDWRKR*M
  2496.                      IWORK = ITAU + M
  2497. *
  2498. *                    Compute A=L*Q
  2499. *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  2500. *
  2501.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2502.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2503. *
  2504. *                    Copy L to WORK(IU), zeroing out below it
  2505. *
  2506.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IU ),
  2507.      $                            LDWRKU )
  2508.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2509.      $                            WORK( IU+LDWRKU ), LDWRKU )
  2510. *
  2511. *                    Generate Q in A
  2512. *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  2513. *
  2514.                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2515.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2516.                      IE = ITAU
  2517.                      ITAUQ = IE + M
  2518.                      ITAUP = ITAUQ + M
  2519.                      IWORK = ITAUP + M
  2520. *
  2521. *                    Bidiagonalize L in WORK(IU), copying result to
  2522. *                    WORK(IR)
  2523. *                    (Workspace: need 2*M*M+4*M,
  2524. *                                prefer 2*M*M+3*M+2*M*NB)
  2525. *
  2526.                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
  2527.      $                            WORK( IE ), WORK( ITAUQ ),
  2528.      $                            WORK( ITAUP ), WORK( IWORK ),
  2529.      $                            LWORK-IWORK+1, IERR )
  2530.                      CALL DLACPY( 'L', M, M, WORK( IU ), LDWRKU,
  2531.      $                            WORK( IR ), LDWRKR )
  2532. *
  2533. *                    Generate right bidiagonalizing vectors in WORK(IU)
  2534. *                    (Workspace: need 2*M*M+4*M-1,
  2535. *                                prefer 2*M*M+3*M+(M-1)*NB)
  2536. *
  2537.                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  2538.      $                            WORK( ITAUP ), WORK( IWORK ),
  2539.      $                            LWORK-IWORK+1, IERR )
  2540. *
  2541. *                    Generate left bidiagonalizing vectors in WORK(IR)
  2542. *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
  2543. *
  2544.                      CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
  2545.      $                            WORK( ITAUQ ), WORK( IWORK ),
  2546.      $                            LWORK-IWORK+1, IERR )
  2547.                      IWORK = IE + M
  2548. *
  2549. *                    Perform bidiagonal QR iteration, computing left
  2550. *                    singular vectors of L in WORK(IR) and computing
  2551. *                    right singular vectors of L in WORK(IU)
  2552. *                    (Workspace: need 2*M*M+BDSPAC)
  2553. *
  2554.                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
  2555.      $                            WORK( IU ), LDWRKU, WORK( IR ),
  2556.      $                            LDWRKR, DUM, 1, WORK( IWORK ), INFO )
  2557. *
  2558. *                    Multiply right singular vectors of L in WORK(IU) by
  2559. *                    Q in A, storing result in VT
  2560. *                    (Workspace: need M*M)
  2561. *
  2562.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
  2563.      $                           LDWRKU, A, LDA, ZERO, VT, LDVT )
  2564. *
  2565. *                    Copy left singular vectors of L to A
  2566. *                    (Workspace: need M*M)
  2567. *
  2568.                      CALL DLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
  2569.      $                            LDA )
  2570. *
  2571.                   ELSE
  2572. *
  2573. *                    Insufficient workspace for a fast algorithm
  2574. *
  2575.                      ITAU = 1
  2576.                      IWORK = ITAU + M
  2577. *
  2578. *                    Compute A=L*Q, copying result to VT
  2579. *                    (Workspace: need 2*M, prefer M+M*NB)
  2580. *
  2581.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2582.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2583.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2584. *
  2585. *                    Generate Q in VT
  2586. *                    (Workspace: need 2*M, prefer M+M*NB)
  2587. *
  2588.                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2589.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2590.                      IE = ITAU
  2591.                      ITAUQ = IE + M
  2592.                      ITAUP = ITAUQ + M
  2593.                      IWORK = ITAUP + M
  2594. *
  2595. *                    Zero out above L in A
  2596. *
  2597.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
  2598.      $                            LDA )
  2599. *
  2600. *                    Bidiagonalize L in A
  2601. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2602. *
  2603.                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
  2604.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2605.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2606. *
  2607. *                    Multiply right vectors bidiagonalizing L by Q in VT
  2608. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  2609. *
  2610.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
  2611.      $                            WORK( ITAUP ), VT, LDVT,
  2612.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2613. *
  2614. *                    Generate left bidiagonalizing vectors of L in A
  2615. *                    (Workspace: need 4*M, prefer 3*M+M*NB)
  2616. *
  2617.                      CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  2618.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2619.                      IWORK = IE + M
  2620. *
  2621. *                    Perform bidiagonal QR iteration, compute left
  2622. *                    singular vectors of A in A and compute right
  2623. *                    singular vectors of A in VT
  2624. *                    (Workspace: need BDSPAC)
  2625. *
  2626.                      CALL DBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
  2627.      $                            LDVT, A, LDA, DUM, 1, WORK( IWORK ),
  2628.      $                            INFO )
  2629. *
  2630.                   END IF
  2631. *
  2632.                ELSE IF( WNTUAS ) THEN
  2633. *
  2634. *                 Path 6t(N much larger than M, JOBU='S' or 'A',
  2635. *                         JOBVT='S')
  2636. *                 M right singular vectors to be computed in VT and
  2637. *                 M left singular vectors to be computed in U
  2638. *
  2639.                   IF( LWORK.GE.M*M+MAX( 4*M, BDSPAC ) ) THEN
  2640. *
  2641. *                    Sufficient workspace for a fast algorithm
  2642. *
  2643.                      IU = 1
  2644.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2645. *
  2646. *                       WORK(IU) is LDA by N
  2647. *
  2648.                         LDWRKU = LDA
  2649.                      ELSE
  2650. *
  2651. *                       WORK(IU) is LDA by M
  2652. *
  2653.                         LDWRKU = M
  2654.                      END IF
  2655.                      ITAU = IU + LDWRKU*M
  2656.                      IWORK = ITAU + M
  2657. *
  2658. *                    Compute A=L*Q
  2659. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2660. *
  2661.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2662.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2663. *
  2664. *                    Copy L to WORK(IU), zeroing out above it
  2665. *
  2666.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IU ),
  2667.      $                            LDWRKU )
  2668.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2669.      $                            WORK( IU+LDWRKU ), LDWRKU )
  2670. *
  2671. *                    Generate Q in A
  2672. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2673. *
  2674.                      CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
  2675.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2676.                      IE = ITAU
  2677.                      ITAUQ = IE + M
  2678.                      ITAUP = ITAUQ + M
  2679.                      IWORK = ITAUP + M
  2680. *
  2681. *                    Bidiagonalize L in WORK(IU), copying result to U
  2682. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  2683. *
  2684.                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
  2685.      $                            WORK( IE ), WORK( ITAUQ ),
  2686.      $                            WORK( ITAUP ), WORK( IWORK ),
  2687.      $                            LWORK-IWORK+1, IERR )
  2688.                      CALL DLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
  2689.      $                            LDU )
  2690. *
  2691. *                    Generate right bidiagonalizing vectors in WORK(IU)
  2692. *                    (Workspace: need M*M+4*M-1,
  2693. *                                prefer M*M+3*M+(M-1)*NB)
  2694. *
  2695.                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  2696.      $                            WORK( ITAUP ), WORK( IWORK ),
  2697.      $                            LWORK-IWORK+1, IERR )
  2698. *
  2699. *                    Generate left bidiagonalizing vectors in U
  2700. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
  2701. *
  2702.                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2703.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2704.                      IWORK = IE + M
  2705. *
  2706. *                    Perform bidiagonal QR iteration, computing left
  2707. *                    singular vectors of L in U and computing right
  2708. *                    singular vectors of L in WORK(IU)
  2709. *                    (Workspace: need M*M+BDSPAC)
  2710. *
  2711.                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
  2712.      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
  2713.      $                            WORK( IWORK ), INFO )
  2714. *
  2715. *                    Multiply right singular vectors of L in WORK(IU) by
  2716. *                    Q in A, storing result in VT
  2717. *                    (Workspace: need M*M)
  2718. *
  2719.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
  2720.      $                           LDWRKU, A, LDA, ZERO, VT, LDVT )
  2721. *
  2722.                   ELSE
  2723. *
  2724. *                    Insufficient workspace for a fast algorithm
  2725. *
  2726.                      ITAU = 1
  2727.                      IWORK = ITAU + M
  2728. *
  2729. *                    Compute A=L*Q, copying result to VT
  2730. *                    (Workspace: need 2*M, prefer M+M*NB)
  2731. *
  2732.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2733.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2734.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2735. *
  2736. *                    Generate Q in VT
  2737. *                    (Workspace: need 2*M, prefer M+M*NB)
  2738. *
  2739.                      CALL DORGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
  2740.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2741. *
  2742. *                    Copy L to U, zeroing out above it
  2743. *
  2744.                      CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
  2745.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
  2746.      $                            LDU )
  2747.                      IE = ITAU
  2748.                      ITAUQ = IE + M
  2749.                      ITAUP = ITAUQ + M
  2750.                      IWORK = ITAUP + M
  2751. *
  2752. *                    Bidiagonalize L in U
  2753. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2754. *
  2755.                      CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
  2756.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2757.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2758. *
  2759. *                    Multiply right bidiagonalizing vectors in U by Q
  2760. *                    in VT
  2761. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  2762. *
  2763.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
  2764.      $                            WORK( ITAUP ), VT, LDVT,
  2765.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2766. *
  2767. *                    Generate left bidiagonalizing vectors in U
  2768. *                    (Workspace: need 4*M, prefer 3*M+M*NB)
  2769. *
  2770.                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  2771.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2772.                      IWORK = IE + M
  2773. *
  2774. *                    Perform bidiagonal QR iteration, computing left
  2775. *                    singular vectors of A in U and computing right
  2776. *                    singular vectors of A in VT
  2777. *                    (Workspace: need BDSPAC)
  2778. *
  2779.                      CALL DBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
  2780.      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
  2781.      $                            INFO )
  2782. *
  2783.                   END IF
  2784. *
  2785.                END IF
  2786. *
  2787.             ELSE IF( WNTVA ) THEN
  2788. *
  2789.                IF( WNTUN ) THEN
  2790. *
  2791. *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
  2792. *                 N right singular vectors to be computed in VT and
  2793. *                 no left singular vectors to be computed
  2794. *
  2795.                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
  2796. *
  2797. *                    Sufficient workspace for a fast algorithm
  2798. *
  2799.                      IR = 1
  2800.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  2801. *
  2802. *                       WORK(IR) is LDA by M
  2803. *
  2804.                         LDWRKR = LDA
  2805.                      ELSE
  2806. *
  2807. *                       WORK(IR) is M by M
  2808. *
  2809.                         LDWRKR = M
  2810.                      END IF
  2811.                      ITAU = IR + LDWRKR*M
  2812.                      IWORK = ITAU + M
  2813. *
  2814. *                    Compute A=L*Q, copying result to VT
  2815. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  2816. *
  2817.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2818.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2819.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2820. *
  2821. *                    Copy L to WORK(IR), zeroing out above it
  2822. *
  2823.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IR ),
  2824.      $                            LDWRKR )
  2825.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2826.      $                            WORK( IR+LDWRKR ), LDWRKR )
  2827. *
  2828. *                    Generate Q in VT
  2829. *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
  2830. *
  2831.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  2832.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2833.                      IE = ITAU
  2834.                      ITAUQ = IE + M
  2835.                      ITAUP = ITAUQ + M
  2836.                      IWORK = ITAUP + M
  2837. *
  2838. *                    Bidiagonalize L in WORK(IR)
  2839. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  2840. *
  2841.                      CALL DGEBRD( M, M, WORK( IR ), LDWRKR, S,
  2842.      $                            WORK( IE ), WORK( ITAUQ ),
  2843.      $                            WORK( ITAUP ), WORK( IWORK ),
  2844.      $                            LWORK-IWORK+1, IERR )
  2845. *
  2846. *                    Generate right bidiagonalizing vectors in WORK(IR)
  2847. *                    (Workspace: need M*M+4*M-1,
  2848. *                                prefer M*M+3*M+(M-1)*NB)
  2849. *
  2850.                      CALL DORGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
  2851.      $                            WORK( ITAUP ), WORK( IWORK ),
  2852.      $                            LWORK-IWORK+1, IERR )
  2853.                      IWORK = IE + M
  2854. *
  2855. *                    Perform bidiagonal QR iteration, computing right
  2856. *                    singular vectors of L in WORK(IR)
  2857. *                    (Workspace: need M*M+BDSPAC)
  2858. *
  2859.                      CALL DBDSQR( 'U', M, M, 0, 0, S, WORK( IE ),
  2860.      $                            WORK( IR ), LDWRKR, DUM, 1, DUM, 1,
  2861.      $                            WORK( IWORK ), INFO )
  2862. *
  2863. *                    Multiply right singular vectors of L in WORK(IR) by
  2864. *                    Q in VT, storing result in A
  2865. *                    (Workspace: need M*M)
  2866. *
  2867.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IR ),
  2868.      $                           LDWRKR, VT, LDVT, ZERO, A, LDA )
  2869. *
  2870. *                    Copy right singular vectors of A from A to VT
  2871. *
  2872.                      CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
  2873. *
  2874.                   ELSE
  2875. *
  2876. *                    Insufficient workspace for a fast algorithm
  2877. *
  2878.                      ITAU = 1
  2879.                      IWORK = ITAU + M
  2880. *
  2881. *                    Compute A=L*Q, copying result to VT
  2882. *                    (Workspace: need 2*M, prefer M+M*NB)
  2883. *
  2884.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2885.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2886.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2887. *
  2888. *                    Generate Q in VT
  2889. *                    (Workspace: need M+N, prefer M+N*NB)
  2890. *
  2891.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  2892.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2893.                      IE = ITAU
  2894.                      ITAUQ = IE + M
  2895.                      ITAUP = ITAUQ + M
  2896.                      IWORK = ITAUP + M
  2897. *
  2898. *                    Zero out above L in A
  2899. *
  2900.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
  2901.      $                            LDA )
  2902. *
  2903. *                    Bidiagonalize L in A
  2904. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  2905. *
  2906.                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
  2907.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  2908.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2909. *
  2910. *                    Multiply right bidiagonalizing vectors in A by Q
  2911. *                    in VT
  2912. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  2913. *
  2914.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
  2915.      $                            WORK( ITAUP ), VT, LDVT,
  2916.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2917.                      IWORK = IE + M
  2918. *
  2919. *                    Perform bidiagonal QR iteration, computing right
  2920. *                    singular vectors of A in VT
  2921. *                    (Workspace: need BDSPAC)
  2922. *
  2923.                      CALL DBDSQR( 'U', M, N, 0, 0, S, WORK( IE ), VT,
  2924.      $                            LDVT, DUM, 1, DUM, 1, WORK( IWORK ),
  2925.      $                            INFO )
  2926. *
  2927.                   END IF
  2928. *
  2929.                ELSE IF( WNTUO ) THEN
  2930. *
  2931. *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
  2932. *                 N right singular vectors to be computed in VT and
  2933. *                 M left singular vectors to be overwritten on A
  2934. *
  2935.                   IF( LWORK.GE.2*M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
  2936. *
  2937. *                    Sufficient workspace for a fast algorithm
  2938. *
  2939.                      IU = 1
  2940.                      IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
  2941. *
  2942. *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
  2943. *
  2944.                         LDWRKU = LDA
  2945.                         IR = IU + LDWRKU*M
  2946.                         LDWRKR = LDA
  2947.                      ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
  2948. *
  2949. *                       WORK(IU) is LDA by M and WORK(IR) is M by M
  2950. *
  2951.                         LDWRKU = LDA
  2952.                         IR = IU + LDWRKU*M
  2953.                         LDWRKR = M
  2954.                      ELSE
  2955. *
  2956. *                       WORK(IU) is M by M and WORK(IR) is M by M
  2957. *
  2958.                         LDWRKU = M
  2959.                         IR = IU + LDWRKU*M
  2960.                         LDWRKR = M
  2961.                      END IF
  2962.                      ITAU = IR + LDWRKR*M
  2963.                      IWORK = ITAU + M
  2964. *
  2965. *                    Compute A=L*Q, copying result to VT
  2966. *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
  2967. *
  2968.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  2969.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2970.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  2971. *
  2972. *                    Generate Q in VT
  2973. *                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
  2974. *
  2975.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  2976.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  2977. *
  2978. *                    Copy L to WORK(IU), zeroing out above it
  2979. *
  2980.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IU ),
  2981.      $                            LDWRKU )
  2982.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  2983.      $                            WORK( IU+LDWRKU ), LDWRKU )
  2984.                      IE = ITAU
  2985.                      ITAUQ = IE + M
  2986.                      ITAUP = ITAUQ + M
  2987.                      IWORK = ITAUP + M
  2988. *
  2989. *                    Bidiagonalize L in WORK(IU), copying result to
  2990. *                    WORK(IR)
  2991. *                    (Workspace: need 2*M*M+4*M,
  2992. *                                prefer 2*M*M+3*M+2*M*NB)
  2993. *
  2994.                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
  2995.      $                            WORK( IE ), WORK( ITAUQ ),
  2996.      $                            WORK( ITAUP ), WORK( IWORK ),
  2997.      $                            LWORK-IWORK+1, IERR )
  2998.                      CALL DLACPY( 'L', M, M, WORK( IU ), LDWRKU,
  2999.      $                            WORK( IR ), LDWRKR )
  3000. *
  3001. *                    Generate right bidiagonalizing vectors in WORK(IU)
  3002. *                    (Workspace: need 2*M*M+4*M-1,
  3003. *                                prefer 2*M*M+3*M+(M-1)*NB)
  3004. *
  3005.                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  3006.      $                            WORK( ITAUP ), WORK( IWORK ),
  3007.      $                            LWORK-IWORK+1, IERR )
  3008. *
  3009. *                    Generate left bidiagonalizing vectors in WORK(IR)
  3010. *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
  3011. *
  3012.                      CALL DORGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
  3013.      $                            WORK( ITAUQ ), WORK( IWORK ),
  3014.      $                            LWORK-IWORK+1, IERR )
  3015.                      IWORK = IE + M
  3016. *
  3017. *                    Perform bidiagonal QR iteration, computing left
  3018. *                    singular vectors of L in WORK(IR) and computing
  3019. *                    right singular vectors of L in WORK(IU)
  3020. *                    (Workspace: need 2*M*M+BDSPAC)
  3021. *
  3022.                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
  3023.      $                            WORK( IU ), LDWRKU, WORK( IR ),
  3024.      $                            LDWRKR, DUM, 1, WORK( IWORK ), INFO )
  3025. *
  3026. *                    Multiply right singular vectors of L in WORK(IU) by
  3027. *                    Q in VT, storing result in A
  3028. *                    (Workspace: need M*M)
  3029. *
  3030.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
  3031.      $                           LDWRKU, VT, LDVT, ZERO, A, LDA )
  3032. *
  3033. *                    Copy right singular vectors of A from A to VT
  3034. *
  3035.                      CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
  3036. *
  3037. *                    Copy left singular vectors of A from WORK(IR) to A
  3038. *
  3039.                      CALL DLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
  3040.      $                            LDA )
  3041. *
  3042.                   ELSE
  3043. *
  3044. *                    Insufficient workspace for a fast algorithm
  3045. *
  3046.                      ITAU = 1
  3047.                      IWORK = ITAU + M
  3048. *
  3049. *                    Compute A=L*Q, copying result to VT
  3050. *                    (Workspace: need 2*M, prefer M+M*NB)
  3051. *
  3052.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  3053.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3054.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3055. *
  3056. *                    Generate Q in VT
  3057. *                    (Workspace: need M+N, prefer M+N*NB)
  3058. *
  3059.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3060.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3061.                      IE = ITAU
  3062.                      ITAUQ = IE + M
  3063.                      ITAUP = ITAUQ + M
  3064.                      IWORK = ITAUP + M
  3065. *
  3066. *                    Zero out above L in A
  3067. *
  3068.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ),
  3069.      $                            LDA )
  3070. *
  3071. *                    Bidiagonalize L in A
  3072. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  3073. *
  3074.                      CALL DGEBRD( M, M, A, LDA, S, WORK( IE ),
  3075.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  3076.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3077. *
  3078. *                    Multiply right bidiagonalizing vectors in A by Q
  3079. *                    in VT
  3080. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  3081. *
  3082.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, A, LDA,
  3083.      $                            WORK( ITAUP ), VT, LDVT,
  3084.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3085. *
  3086. *                    Generate left bidiagonalizing vectors in A
  3087. *                    (Workspace: need 4*M, prefer 3*M+M*NB)
  3088. *
  3089.                      CALL DORGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
  3090.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3091.                      IWORK = IE + M
  3092. *
  3093. *                    Perform bidiagonal QR iteration, computing left
  3094. *                    singular vectors of A in A and computing right
  3095. *                    singular vectors of A in VT
  3096. *                    (Workspace: need BDSPAC)
  3097. *
  3098.                      CALL DBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
  3099.      $                            LDVT, A, LDA, DUM, 1, WORK( IWORK ),
  3100.      $                            INFO )
  3101. *
  3102.                   END IF
  3103. *
  3104.                ELSE IF( WNTUAS ) THEN
  3105. *
  3106. *                 Path 9t(N much larger than M, JOBU='S' or 'A',
  3107. *                         JOBVT='A')
  3108. *                 N right singular vectors to be computed in VT and
  3109. *                 M left singular vectors to be computed in U
  3110. *
  3111.                   IF( LWORK.GE.M*M+MAX( N+M, 4*M, BDSPAC ) ) THEN
  3112. *
  3113. *                    Sufficient workspace for a fast algorithm
  3114. *
  3115.                      IU = 1
  3116.                      IF( LWORK.GE.WRKBL+LDA*M ) THEN
  3117. *
  3118. *                       WORK(IU) is LDA by M
  3119. *
  3120.                         LDWRKU = LDA
  3121.                      ELSE
  3122. *
  3123. *                       WORK(IU) is M by M
  3124. *
  3125.                         LDWRKU = M
  3126.                      END IF
  3127.                      ITAU = IU + LDWRKU*M
  3128.                      IWORK = ITAU + M
  3129. *
  3130. *                    Compute A=L*Q, copying result to VT
  3131. *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
  3132. *
  3133.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  3134.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3135.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3136. *
  3137. *                    Generate Q in VT
  3138. *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
  3139. *
  3140.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3141.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3142. *
  3143. *                    Copy L to WORK(IU), zeroing out above it
  3144. *
  3145.                      CALL DLACPY( 'L', M, M, A, LDA, WORK( IU ),
  3146.      $                            LDWRKU )
  3147.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
  3148.      $                            WORK( IU+LDWRKU ), LDWRKU )
  3149.                      IE = ITAU
  3150.                      ITAUQ = IE + M
  3151.                      ITAUP = ITAUQ + M
  3152.                      IWORK = ITAUP + M
  3153. *
  3154. *                    Bidiagonalize L in WORK(IU), copying result to U
  3155. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
  3156. *
  3157.                      CALL DGEBRD( M, M, WORK( IU ), LDWRKU, S,
  3158.      $                            WORK( IE ), WORK( ITAUQ ),
  3159.      $                            WORK( ITAUP ), WORK( IWORK ),
  3160.      $                            LWORK-IWORK+1, IERR )
  3161.                      CALL DLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
  3162.      $                            LDU )
  3163. *
  3164. *                    Generate right bidiagonalizing vectors in WORK(IU)
  3165. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
  3166. *
  3167.                      CALL DORGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
  3168.      $                            WORK( ITAUP ), WORK( IWORK ),
  3169.      $                            LWORK-IWORK+1, IERR )
  3170. *
  3171. *                    Generate left bidiagonalizing vectors in U
  3172. *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
  3173. *
  3174.                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  3175.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3176.                      IWORK = IE + M
  3177. *
  3178. *                    Perform bidiagonal QR iteration, computing left
  3179. *                    singular vectors of L in U and computing right
  3180. *                    singular vectors of L in WORK(IU)
  3181. *                    (Workspace: need M*M+BDSPAC)
  3182. *
  3183.                      CALL DBDSQR( 'U', M, M, M, 0, S, WORK( IE ),
  3184.      $                            WORK( IU ), LDWRKU, U, LDU, DUM, 1,
  3185.      $                            WORK( IWORK ), INFO )
  3186. *
  3187. *                    Multiply right singular vectors of L in WORK(IU) by
  3188. *                    Q in VT, storing result in A
  3189. *                    (Workspace: need M*M)
  3190. *
  3191.                      CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IU ),
  3192.      $                           LDWRKU, VT, LDVT, ZERO, A, LDA )
  3193. *
  3194. *                    Copy right singular vectors of A from A to VT
  3195. *
  3196.                      CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
  3197. *
  3198.                   ELSE
  3199. *
  3200. *                    Insufficient workspace for a fast algorithm
  3201. *
  3202.                      ITAU = 1
  3203.                      IWORK = ITAU + M
  3204. *
  3205. *                    Compute A=L*Q, copying result to VT
  3206. *                    (Workspace: need 2*M, prefer M+M*NB)
  3207. *
  3208.                      CALL DGELQF( M, N, A, LDA, WORK( ITAU ),
  3209.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3210.                      CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3211. *
  3212. *                    Generate Q in VT
  3213. *                    (Workspace: need M+N, prefer M+N*NB)
  3214. *
  3215.                      CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
  3216.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3217. *
  3218. *                    Copy L to U, zeroing out above it
  3219. *
  3220.                      CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
  3221.                      CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, U( 1, 2 ),
  3222.      $                            LDU )
  3223.                      IE = ITAU
  3224.                      ITAUQ = IE + M
  3225.                      ITAUP = ITAUQ + M
  3226.                      IWORK = ITAUP + M
  3227. *
  3228. *                    Bidiagonalize L in U
  3229. *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
  3230. *
  3231.                      CALL DGEBRD( M, M, U, LDU, S, WORK( IE ),
  3232.      $                            WORK( ITAUQ ), WORK( ITAUP ),
  3233.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3234. *
  3235. *                    Multiply right bidiagonalizing vectors in U by Q
  3236. *                    in VT
  3237. *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
  3238. *
  3239.                      CALL DORMBR( 'P', 'L', 'T', M, N, M, U, LDU,
  3240.      $                            WORK( ITAUP ), VT, LDVT,
  3241.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3242. *
  3243. *                    Generate left bidiagonalizing vectors in U
  3244. *                    (Workspace: need 4*M, prefer 3*M+M*NB)
  3245. *
  3246.                      CALL DORGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
  3247.      $                            WORK( IWORK ), LWORK-IWORK+1, IERR )
  3248.                      IWORK = IE + M
  3249. *
  3250. *                    Perform bidiagonal QR iteration, computing left
  3251. *                    singular vectors of A in U and computing right
  3252. *                    singular vectors of A in VT
  3253. *                    (Workspace: need BDSPAC)
  3254. *
  3255.                      CALL DBDSQR( 'U', M, N, M, 0, S, WORK( IE ), VT,
  3256.      $                            LDVT, U, LDU, DUM, 1, WORK( IWORK ),
  3257.      $                            INFO )
  3258. *
  3259.                   END IF
  3260. *
  3261.                END IF
  3262. *
  3263.             END IF
  3264. *
  3265.          ELSE
  3266. *
  3267. *           N .LT. MNTHR
  3268. *
  3269. *           Path 10t(N greater than M, but not much larger)
  3270. *           Reduce to bidiagonal form without LQ decomposition
  3271. *
  3272.             IE = 1
  3273.             ITAUQ = IE + M
  3274.             ITAUP = ITAUQ + M
  3275.             IWORK = ITAUP + M
  3276. *
  3277. *           Bidiagonalize A
  3278. *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
  3279. *
  3280.             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
  3281.      $                   WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
  3282.      $                   IERR )
  3283.             IF( WNTUAS ) THEN
  3284. *
  3285. *              If left singular vectors desired in U, copy result to U
  3286. *              and generate left bidiagonalizing vectors in U
  3287. *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
  3288. *
  3289.                CALL DLACPY( 'L', M, M, A, LDA, U, LDU )
  3290.                CALL DORGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
  3291.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3292.             END IF
  3293.             IF( WNTVAS ) THEN
  3294. *
  3295. *              If right singular vectors desired in VT, copy result to
  3296. *              VT and generate right bidiagonalizing vectors in VT
  3297. *              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
  3298. *
  3299.                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
  3300.                IF( WNTVA )
  3301.      $            NRVT = N
  3302.                IF( WNTVS )
  3303.      $            NRVT = M
  3304.                CALL DORGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
  3305.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3306.             END IF
  3307.             IF( WNTUO ) THEN
  3308. *
  3309. *              If left singular vectors desired in A, generate left
  3310. *              bidiagonalizing vectors in A
  3311. *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
  3312. *
  3313.                CALL DORGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
  3314.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3315.             END IF
  3316.             IF( WNTVO ) THEN
  3317. *
  3318. *              If right singular vectors desired in A, generate right
  3319. *              bidiagonalizing vectors in A
  3320. *              (Workspace: need 4*M, prefer 3*M+M*NB)
  3321. *
  3322.                CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
  3323.      $                      WORK( IWORK ), LWORK-IWORK+1, IERR )
  3324.             END IF
  3325.             IWORK = IE + M
  3326.             IF( WNTUAS .OR. WNTUO )
  3327.      $         NRU = M
  3328.             IF( WNTUN )
  3329.      $         NRU = 0
  3330.             IF( WNTVAS .OR. WNTVO )
  3331.      $         NCVT = N
  3332.             IF( WNTVN )
  3333.      $         NCVT = 0
  3334.             IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
  3335. *
  3336. *              Perform bidiagonal QR iteration, if desired, computing
  3337. *              left singular vectors in U and computing right singular
  3338. *              vectors in VT
  3339. *              (Workspace: need BDSPAC)
  3340. *
  3341.                CALL DBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), VT,
  3342.      $                      LDVT, U, LDU, DUM, 1, WORK( IWORK ), INFO )
  3343.             ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
  3344. *
  3345. *              Perform bidiagonal QR iteration, if desired, computing
  3346. *              left singular vectors in U and computing right singular
  3347. *              vectors in A
  3348. *              (Workspace: need BDSPAC)
  3349. *
  3350.                CALL DBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), A, LDA,
  3351.      $                      U, LDU, DUM, 1, WORK( IWORK ), INFO )
  3352.             ELSE
  3353. *
  3354. *              Perform bidiagonal QR iteration, if desired, computing
  3355. *              left singular vectors in A and computing right singular
  3356. *              vectors in VT
  3357. *              (Workspace: need BDSPAC)
  3358. *
  3359.                CALL DBDSQR( 'L', M, NCVT, NRU, 0, S, WORK( IE ), VT,
  3360.      $                      LDVT, A, LDA, DUM, 1, WORK( IWORK ), INFO )
  3361.             END IF
  3362. *
  3363.          END IF
  3364. *
  3365.       END IF
  3366. *
  3367. *     If DBDSQR failed to converge, copy unconverged superdiagonals
  3368. *     to WORK( 2:MINMN )
  3369. *
  3370.       IF( INFO.NE.0 ) THEN
  3371.          IF( IE.GT.2 ) THEN
  3372.             DO 50 I = 1, MINMN - 1
  3373.                WORK( I+1 ) = WORK( I+IE-1 )
  3374.    50       CONTINUE
  3375.          END IF
  3376.          IF( IE.LT.2 ) THEN
  3377.             DO 60 I = MINMN - 1, 1, -1
  3378.                WORK( I+1 ) = WORK( I+IE-1 )
  3379.    60       CONTINUE
  3380.          END IF
  3381.       END IF
  3382. *
  3383. *     Undo scaling if necessary
  3384. *
  3385.       IF( ISCL.EQ.1 ) THEN
  3386.          IF( ANRM.GT.BIGNUM )
  3387.      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
  3388.      $                   IERR )
  3389.          IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
  3390.      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, WORK( 2 ),
  3391.      $                   MINMN, IERR )
  3392.          IF( ANRM.LT.SMLNUM )
  3393.      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
  3394.      $                   IERR )
  3395.          IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
  3396.      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, WORK( 2 ),
  3397.      $                   MINMN, IERR )
  3398.       END IF
  3399. *
  3400. *     Return optimal workspace in WORK(1)
  3401. *
  3402.       WORK( 1 ) = MAXWRK
  3403. *
  3404.       RETURN
  3405. *
  3406. *     End of DGESVD
  3407. *
  3408.       END
  3409.